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A  FLIGHT  DYNAMIC  SIMULATION  PROGRAM 
IN  AIR-PATH  AXES  USING  A.C.S.L. 

by 

P.W.  G  IBSENS 

SUMMARY 

The  six  degrees  of  freedom  dynamic  equations  of  motion  have  been 
programmed  in  the  Advanced  Continuous  Simulation  Language  (ACSL)  for  use 
in  aircraft  simulations  at  ARL.  Air-path  axes  were  chosen  for  the  integration 
of  the  force  equations,  and  body  axes  for  the  integration  of  the  moment 
equations.  The  use  of  quaternions  for  the  determination  of  the  direction 
cosines  has  been  described. 
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NOTATION 


Linear  accelerations  in  body  axes  (g’s) 

Normal  acceleration  (g's) 

Inertia  constants 

Non-dimensional  aerodynamic  coefficients 

Equilibrium  drag  coefficient 

Derivative  of  drag  coefficient  with  respect  to  velocity 

Derivative  of  drag  coefficient  with  respect  to  angle  of 

attack 

Equilibrium  lift  coefficient 

Derivative  of  lift  coefficient  with  respect  to  pitch  rate 
Derivative  of  lift  coefficient  with  respect  to  velocity 

Derivative  of  lift  coefficient  with  respect  to  angle  of 

attack 

Derivative  of  lift  coefficient  with  respect  to  rate  of 
change  of  angle  of  attack 

Derivative  of  moment  coefficient  with  respect  to  pitch  rate 
Derivative  of  moment  coefficient  with  respect  to  velocity 
Derivative  of  moment  coefficient  with  respect  to  angle  of 
attack 

Derivative  of  moment  coefficient  with  respect  to  rate  of 
change  of  angle  of  attack 
Equilibrium  thrust  coefficient 

Derivative  of  thrust  coefficient  with  respect  to  velocity 
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NOTATION  (Contd.) 


GX* 

FGY’ 

fgz 

Gravitational  forces 

(N) 

XB* 

fyb» 

fzb 

Total  forces  in  body  axes 

(N) 

XS* 

fys* 

Fzs 

Total  forces  in  stability  axes 

(N) 

'xw 

fyw* 

Fzw 

Total  forces  in  air-path  axes 

(N) 

Gravitational  constant 

(ms' 

XX* 

IYY» 

*zz*  Ixz 

Moments  and  product  of  inertia 

(Kg 

L,  M,  N  Total  applied  moments  in  body  axes  (Nm) 
lA*  mA'  na  Applied  aerodynamic  moments  in  stability 

axes  (Nm) 
Lp,  Mp,  Np  Applied  propulsive  moments  in  body  axes(Nm) 
Li ,  (i=1,3)  Direction  cosines 


P,  Q,  R 

PS»  RS 
RANGE 

UB*  VB’  ^B 
VE  *  VEK 
VGR 

VNE*  VEE  *  VDE 


Aircraft  mass  (Kg) 

Angular  velocity  components  about  body 
axes  (rad.  s-1) 

Pitch  and  yaw  rates  about  stability  axes(rad.s~1 ) 
Aircraft  range  (m) 

Body  axes  velocity  components  (ms  1 ) 


Equivalent  airspeed  (ms  ,  Knots) 

Aircraft  ground  speed  (ms-1) 

Components  of  velocity  relative  to  the  earth(ms-1) 
True  airspeed  (ms  ') 

Applied  aeroaynamic  forces  in  stability  axes (N ) 
Applied  propulsive  forces  in  body  axes  (N) 

Body  axes  reference  frame 
Earth  axes  reference  frame 

Stability  axes  reference  frame 
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NOTATION  (Contd.) 


Xw,  ¥w,  Zw  Air-path  axes  reference  frame 

a  Angle  of  attack  (rad) 

3  Angle  of  sideslip  (rad) 

Y  Flight  path  angle  (rad) 

0,  <j>,  i|>  Euler  angles  of  pitch,  roll  (bank)  and  yaw 

(heading)  (rad) 

X  Angle  of  climb  (rad) 

t  (i=0,3)  Quaternion  components 

X  Angle  of  track,  east  of  north  (rad) 


Subscripts 

A 

B 

E 

P 

S 

W 


Aerodynamic  contribution 
Body  axes 
Earth  axes 

Propulsive  contribution 
Stability  axes 
Air-path  axes 


A  dot  over  a  variable  denotes  the  first  derivative  with  respect  to  time 
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1 .  INTRODUCTION 

The  aircraft  behavioural  Studies  -  Fixed  Wing  (ABS  -  FW)  group 
of  Aerodynamics  Division  at  the  Aeronautical  Research  Laboratories  (ARL) 
is  concerned  with  the  flight  dynamic  behaviour  of  fixed  wing  aircraft,  and 
has  the  responsibility  for  developing  flight  dynamic  computer  models  of 
the  advanced  nigh  performance  aircraft  operated  by  the  RAAF . 

This  memorandum  documents  the  basic  six  degrees  of  freedom 
dynamic  equations  of  motion  incorporated  in  the  associated  simulation 
program  SDOFAP  which  was  written  using  Advanced  Continuous  Simulation 
Language  (ACSL).  The  program  has  been  developed  for  general  use  by  the 
ABS-FW  group  for  aircraft  simulation  studies.  The  program  presented  in 
Ref.  (1),  which  was  originally  written  using  Earth  axes,  has  been  modified 
to  use  Air-path  axes  for  integration  of  the  force  equations  to  allow  the 
linear  analysis  capabilities  of  ACSL  to  be  utilised  more  conveniently. 

Section  2  of  the  memorandum  give3  a  description  of  the  program 
and  its  structure,  while  section  3  deals  with  the  axes  systems,  their 
selection  and  transformation ,  the  use  of  quaternions  in  determining 
aircraft  attitude,  and  the  six  degrees  of  freedom  equations  of  motion. 
Section  4  contains  an  example  application  of  the  program  which 
demonstrates  the  analytical  capaoilities  of  the  ACSL  language  and 
illustrates  various  alternative  presentations  of  the  output  data. 


2.  ACSL  AIRCRAFT  FLIGHT  DYNAMIC  SIMULATION  PROGRAM  -  SDOFAP 

Tne  six  degree  of  freedom  aircraft  simulation  program  SDOFAP  has 
been  written  in  air-path  axes  and  utilises  tne  analytical  features  of  the 
ACSL  language.  A  description  of  ACSL  and  its  utilisation  is  given  in  Kef. 


The  ACSL  program  structure  contains  three  primary  regions,  each 
dealing  with  specific  segments  of  program  execution;  these  are  the 
INITIAL,  DYNAMIC  and  DERIVATIVE  sections.  Fig.  1  illustrates  the  program 
flow. 

The  INITIAL  section  serves  the  following  purposes; 

1  .  Aircraft  configuration  and  aerodynamic  data  are  read  in  at 
subroutine  level  from  an  input  file.  This  file  is  prepared 
by  the  setup  program  FTCHOO. 

2.  The  initialisation  of  state  and  control  variables  and 
specification  of  atmospheric  data  and  other  constants  is 
performed . 

3.  The  setting  and  updating  of  prescribed  control  variables  and 
runtime  parameters  for  particular  simulation  runs  is 
performed . 

Time  history  calculation  is  carried  out  within  the  DYNAMIC 
section.  This  portion  of  the  program  primarily  administers  the  trimming 
procedures,  prepares  output  variables  and  manages  the  generation  and 
logging  of  time  dependent  results. 

Tne  DERIVATIVE  section  contains  tne  detuls  of  the  six  degrees 
of  freedom  flight  model. 


Tr imming 

of  the  aircraft 

equations  is  peri'or 
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supplied  trimming 

subroutine.  The 

routine  PJwlT  used 

in  the 

example 

application  calls 

the  subroutine 
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common 

statements,  EVAL  must  be  appended  to  the  ACSL  program  and  variable  values 
made  available  by  use  of  the  ACSL  inclusion  character  as  described  in 


Ref.  (2). 


To  provide  the  trim  routine  with  initial  state  and  control 
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variables,  an  approximation  is  made  Defore  trimming.  When  trimming  is 
complete,  control  is  transferred  back  to  the  INITIAL  section. 

A  description  of  the  six  Gegree  of  freedom  equations  is  given  in 
Section  3»  and  complete  program  listings  are  presented  in  Appendix  1. 


3.  SDOFAP  SIMULATION  MODEL 
3.1  Definition  of  Axes  Systems 


m  - 


All  axes  systems  are  assumed  to  be  orthogonal,  right-handed 


.  • 


triads,  and  are  shown  in  figure  2. 


(i)  Earth  Axes  (XE,YE,Z£) 

The  origin  is  at  a  point  fixed  on  the  earth's  surface,  typically 
at  the  runway  threshold  and  on  the  centreline.  The  x-axis  points  North, 
the  Y-axis  East,  and  the  Z-axis  'down'  toward  the  centre  of  the  eartn.  It 
is  assumed  that  the  eartn  is  flat  anc  non-rotatinb,  such  that  the  earth 
axes  are  regarded  as  an  inertial  frame. 


(ii)  Body  Axes  (XB,YB,Zb) 

douy  axes  are  fixed  on  the  aircraft  witn  the  origin  loc  it--.,  at 
tne  aircraft  centre  of  gravity.  Tne  aircraft  is  assumed  to  oe  rig; i,  w .  t  n 
the  X-axis  parallel  to  the  norinonta,.  fuselage  reference  line  ana  p,inting 
'forward',  tne  Y-axis  pointing  to  starboard  trignt),  (,;io  the  g-  i  . 

'  downward ' . 
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(iii)  Stability  Axes  (X^.Yg.Zg) 

Stability  axes  are  a  special  set  of  body  axes  used  primarily  in 
the  study  of  small  disturbances  from  a  steady  reference  flight 
condition.  Aerodynamic  data  are  frequently  presented  in  stability  axes. 
These  axes  are  displaced  from  the  body  frame  by  the  angle  of  attack,  a, 
such  that  the  X-axis  in  the  steady-state  is  aligned  with  the  projection  of 
the  relative  wind  vector  on  the  aircraft  plane  of  symmetry;  the  Y-axis 
points  to  starboard,  and  the  Z-axis  'downward'. 


(iv)  Air-Path  Axes  (XW,YW,ZW) 

Air-path  axes  differ  from  oody  axes  by  the  angle  of 
attack,  a,  and  the  angle  of  sideslip,  6.  The  transformation  from  body  to 
air-path  axes  is  accomplished  as  shown  in  figure  2,  by  first  pitching 
through  -a,  to  coincide  with  the  stability  axes,  and  then  yawing 
through  0.  The  origin  is  located  at  the  aircraft  centre  of  gravity  and 
the  X-axis  is  aligned  with  the  relative  wind  vector. 


Note:  When  tne  wind  velocity  components  are  oero,  the  air-patn  axes 

coincide  with  the  flight-path  axes  as  defined  by  Fogarty  and  Howe  (H). 
Etkin  (3)  refers  to  tne  air-path  axes  as  the  air-trajectory  reference 
frame  (or  wind  axes). 


3.2  Selection  of  Axes 

Tne  various  options  for  selection  of  axes  systems  are  discussed 
in  Ref.  (3)  and  Ref.  1*0.  with  the  development  of  nigh  speed  digital 
computers  with  large  data  storage,  t.-.c-  selection  of  axes  systems  nas 
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become  less  critical. 
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(i)  Force  Equations. 

In  recent  years,  the  use  of  body  axes  for  the  computation  and 
integration  of  force  equations  has  become  less  popular.  Ref.  (1),  from 
which  this  program  was  developed,  employs  earth  axes  to  suit  the 
particular  application  in  that  report.  However,  air- path  axes  have  been 
selected  here  because  when  used  in  conjunction  with  the  ACSL  linear 
analysis  procedures  they  yield  parameters  which  can  be  used  directly  for 
simulation  validation  and  assessment. 

(ii)  Moment  Equations 

The  body  axes  system  is  the  natural  choice  for  the  solution  of 
the  rotational  equations  of  motion  because  of  the  important  advantage  of 
constant  moments  of  inertia  when  calculating  the  moments  and  angular 


motion  of  the 

aircraft . 

3-3  Aircraft 

Attitude 

Determination 

The 

attitude 

of  an 

aircraft  is  defined 

in  terms  of  the 

tradi tional 

Euler 

angles,  ip 

(heading  angle),  e 

(pitch  attitude), 

and  <J>  (roll , 

or  'bank' 

angle ) . 

In  order  to  avoid  the 

problems  associated 

with  the  singularity  in  the  Euler  'rate'  equations,  wmch  occurs 
when  e  =±y0°,  quaternion  components  ( b)  or  direction  cosines  may  be  used 
in  the  integration  step. 

Quaternion  components  were  chosen  for  the  following  reasons: 

U)  their  time  derivates  are  always  finite  and  continuous, 
whereas  those  of  tne  Euler  anbles  possess  singularities; 

(ii)  the  computations  remain  accurate  as  0  ipproacnes  jO- ; 


(iii ) 


it  is  a  four  parameter  method  consisting  of  four 


integrations  with  one  constraint  equation,  whereas  direction 
cosines,  in  principle,  require  nine  integrations  and  six 
constraint  equations. 

The  quaternion  components  are  expressible  in  terms  of  Euler 
angles  as  follows: 

t q  =  cos<{>/2  cos6/2  cos’ii/2  +  siniji/2  sin9/2  simji/2 

t  =  sim p/2  cose/2  costii/2  -  cosij>/2  sine/2  simji/2 

=  coS(j)/2  sine/2  cosi|;/2  +  sin<j>/2  cos 0/2  sirup/2 

=  cos<p/2  cos 0/2  simii/2  -  sinij>/2  sin0/2  cosiji/2 


The  quaternion  component  time  derivatives  are  given  by, 


1/2 

(Pt1  ♦ 

Qj2  ♦ 

Rt^) 

1/2 

(Pto  “ 

QT3  + 

Rt  2} 

(2) 

1/2 

(Pt3  + 

Qt0  " 

Rt  ^ ) 

1/2 

( Px  ^  “ 

Qx,  - 

Rtc) 

where  P,Q,R  are  angular  velocity  components  about  body  axes,  ana 


(3) 


Failure  to  normalize  the  quaternion  components  at  each  iteration  c 
result  in  the  integration  becoiiiinb  unstable. 


►  v  v  vv  v  v>  wwv; 


Euler  angles  may  be  derived  from  the  quaternion  components  or 
from  direction  cosines  by  using  the  following  relationships, 


XT  ♦  T  X 
-1  0  1 
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0  =  tan 


_ t0t2  ~  tit3 _ 

2  2^  2  -  2  T/2 

[(t0  +t,  -  1/2)  +(V2+t0T3)  ] 
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The  initial  Euler  angles  are  used  to  determine  the  initial 
quaternion  components,  which  are  in  turn  used  to  calculate  the  direction 
cosine  parameters  for  use  in  axes  transformation  computat ions .  The 
quaternion  components  are  updated  at  each  iteration,  using  equation  (2), 
such  that  the  direction  cosines  are  recalculated  for  use  in  the  equations 
of  motion,  while  the  Euler  angles  are  calculated  as  output  data  only. 


j.4  Axe3  Transformation 

Trarusf ormat ion  of  u  set  of  variables  from  body  axes  t  j  eurfn 
axes  (or  vice-versa)  is  conveniently  achieved  by  use  of  direction  cosines 
(/),  which  are  obt  airieu  in  terms  of  tne  yuatern.ou  component s  by  the 
following  relit ionsni ps , 
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3.5  Equations  of  Motion 

Figure  3  is  a  summary  of  the  overall  six  degrees  of  freedom 
dynamic  equations  for  the  case  of  a  flat  earth. 

(i)  Force  Equations 

The  aerodynamic  force  components  are  frequently  computed  along 
stability  axes,  and  the  propulsive  force  components  are  usually  supplied 
in  body  axes.  Resolution  along  stability  axes  is  given  by: 


FXS  -  XP  C03a  +  Zpsina  +  XA  +  FGX 
fys  -  YP  +  ya  +  fgy 


F..„  =  cosa  -  X„  sina  +  Z.  +  F„„ 
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where  FQX,  FGy»  Fqz  are  the  gravitational  forces  resolved  into  stability 
axes  through  use  of  the  direction  cosines. 

F_v  =  L_,mg  cosa  +  N0mg  sina 
uX  i  3 

Fqy  =  M^mg  (9) 

Fqz  =  -L^mg  sina  +  N^mg  cosa 
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'g'  is  assumed  to  be  constant  such  that  the  calculated  altitude  in  figure 
3  is  the  geopotential  height,  as  used  in  standard  atmosphere  calculations. 
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The  total  force  components  in  air-path  axes  are  obtained  by 
transformation  of  the  forces  in  stability  axes  through  the  sideslip 
angle  6- 


Fxw  3  Fxs  cose  +  fys  sin6 
fyw  “  “  Fxs  slnB  *  fys  cose 
Fzw  =  Fzs 


The  state  variaole  derivatives  are  obtained  from  the  dynamic  equations; 


u  cy  cose  -  PgSing  +  Fzw/mVT)/cosB 
g  =  F„  /mV  -  R 

Hi  1  O 

V„  =  F. ,,,/m 
T  aW 


where  m  is  the  aircraft  mass,  and  virtual  mass  effects  are  ignored.  P^ 
and  Rc  are  the  angular  rates  about  tne  X  and  Z  axes  respectively,  in 
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staoility  axes. 


P  »  P  cosa  +  R  sina 

O 

Rs  »  -  P  sina  +  R  cos a 


The  velocity  components  of  the  aircraft  relative  to  the  air  may 
be  obtained  in  body  axes,  if  required,  from  the  state  variables,  after 
integration  of  equation  (11). 


UB  «  VT  cos  a  cos  6 

VB  *  VT  sin  6 

Wa  =  V_  sin  a  cos  6 

D  1 


Wind  components  are  introduced  to  give  the  components  of  aircraft  velocity 
relative  to  the  earth. 
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wnere  are  the  wind  component  s  North,  East  and  Down 

respectively  with  respect  to  the  earth. 

Flight,  path  parameters  are  deriveu  directly  from  the  earth  axes 
velocity  vector  components  using  equations  (1b)  to  MY). 
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Ground  speed 


2  2  V 

V  =  (V^  +  V  )  2 

GR  v  NE  EE' 


Climb  angle , 


»  •  T“  WDe/vgr] 


Angle  of  Track, 


X  =Tan  ^Vee/Vne^ 


also,  equation  (18)  gives  the  flight  path  angle  a. 


7  =  6  -  a 


The  positional  coordinates  of  the  aircrafts  centre  of  gravity  are  deduced 
by  integrating  the  earth  axes  velocity  vector  components: 


XE  =  VNE 
*E  =  VEE 
^E  =  "VDE 


to  give  distance  North,  distance  East  and  Altitude  respectively,  where 


altitude  (ALT) 
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The  additional  parameter,  the  Range  of  the  aircraft  C.G.  from 
the  runway  threshold  is  calculated  using  the  definition: 


RANGE  =  (X£2  +  fE2)^ 


(22) 


The  linear  accelerations  sensed  by  accelerometers  mounted  at  the 
c.g.  and  aligned  along  the  body  axes  are  computed  from  the  applied 
aerodynamic  and  propulsion  forces  as  follows: 


fxb 

"  XA 

cosa 

F„ 

=  Y 

+  Y 

YB 

A 

P 

**3 

CD 

II 

X 

> 

sina 

(23) 


The  linear  accelerations  are  then  given  by 


FXB/mg 

Fyy/mg 

FZB/mg 


(2iJ) 


In  aircraft  operations  reference  is  more  commonly  maae  to  the 
normal  acceleration  or  'load  factor'  which  is  defined  as: 


a  =  -a 
n  ; 


(2b) 


(ii)  Moment  Equations 

Tne  total  moments  acting  on  an  aircraft  consist  of  aerodyna::u 
and  powerpiant  components.  Tne  po  *  o  r  p i u :  1 1  components  whien  mclum 
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gyroscopic  moments  due  to  powerplant  rotors,  and  thrust  alignment  moments, 
are  normally  given  in  body  axes.  The  aerodynamic  moments  are,  like  the 
aerodynamic  forces,  frequently  given  in  stability  axes  (3).  If  the 
aerodynamic  moments  are  given  in  body  axes,  the  simpl if ication  of  equation 
(26)  is  obvious  (i.e.  set  a=0). 

L  =  L.  cosa  -  N,  sina  +  L_ 

A  A  P 

M  =  M  +  Mn  (26) 

A  P 

N  =  sina  +  cosa  +  Np 

If  it  is  assumed  that  the  aircraft  has  a  piane  of  symmetry,  such 
that  the  products  of  inertia  IyZ  and  IXY  are  zero,  then  the  body  axes 
angular  accelerations  can  be  calculated  by  using  equation  (27). 
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The  constants  Cg  to  Cg  are  evaluated  during  initialization. 

The  aircraft  angular  velocity  components  may  then  be  obtained  by 
integrating  equation  (27). 


it.  EXAMPLE  APPLICATION:  SIMULATION  OF  THE  LONGITUDINAL  MOTION 

OF  A  LIGHT  AIRCRAFT 


4.1  Background. 

In  Ref.  (6),  the  effects  of  power  on  the  longitudinal 
aerodynamic  characteristics  of  a  single  engined,  propeller  driven 
aeroplane  were  investigated.  The  simulation  developed  in  that  report  was 
written  in  FORTRAN  66  code  for  the  DEC  system  10  computer.  The 
subroutines  employed  a  model  of  propulsion  effects  on  the  aerodynamics  of 
a  competition  aerobatic  aircraft.  These  subroutines  have  been  updated  to 
FORTRAN  77  level  and  transferred  to  an  ELXSI  6400  computer  for  use  in  the 
SDOFAP  model. 

Appendix  2  presents  a  listing  of  the  simulation  program  using 
the  master  program  SDOFAP. ACSL,  plus  the  associated  subroutines.  Since 
the  aerodynamic  and  propulsive  forces  are  inter-related  in  the  power 
effects  modelling,  tney  have  oeen  supplied  to  the  ACSL  program  already 
combined  in  stability  axes,  through  the  subroutine  AERO,  and  the 
subroutine  PROP  which  normally  supplies  the  propulsive  forces  was 
bypassed . 


The  calling  sequence  of  tne  power  effects  subroutines  is 
illustrated  in  Fit,.  4  togefner  witn  crief  explanations  of  their  individual 


purposes . 


Preparation  of  the  input  data  for  this  example  is  demonstrated 
in  Appendix  3  together  with  a  listing  of  the  data  preparation  program 
FTCHOO.  Some  of  the  variables  listed  in  the  data  tables  relate  to  the 
original  power  effects  simulation  program  of  Ref.  (6)  and  are  not  used  in 
the  ACSL  program.  It  is  recommended  however  that  this  format  be  employed 
for  data  entry. 

4.2  Time  Histories 

Time  histories  of  the  longitudinal  dynamic  behaviour  of  the 
light  aircraft  have  been  generated  in  Appendix  4.  Input  data  and  command 
files  have  been  included  to  demonstrate  the  use  of  ACSL  run-time  commands. 

Both  the  short  period  and  phugoid  modes  have  been  analysed  to 
illustrate  the  various  plotting  capabilities  of  ACSL. 

4.3  Eigen  Analysis 

The  eigen  analysis  of  Appendix  b  gives  results  for  the  same 
flight  case  as  in  section  4.2. 

By  eliminating  the  lateral  variables  using  the  FREEZE  facility, 
results  for  the  longitudinal  modes  only  are  ootained.  The  first 
eigenvalue  and  its  accompanying  eigenvector  are  associated  with  the 
quaternions,  the  second  and  third  with  the  phugoid  mode  and  the  fourth  and 
fifth  with  tiie  short  period  mode. 

4.4  Jacobian  Analysis 

The  purpose  of  this  extension  to  tne  program  is  to  extract  tne 
non-dimensional  aerodynamic  derivatives  from  tne  elements  of  the  non- 
dimensional  Jacobian  matrix. 


1  VI-.7  V" 


(  16) 


The  ANALYZ  'JACOB'  run-time  command  causes  ACSL  to  calculate  a 

Jacobian  matrix  around  the  current  trim  point  in  state  space.  In  order  to 

gain  access  to  this  matrix  it  is  necessary  to  call  the  subroutine  INTERM 

from  the  ACSL  library  subroutine  ZZEIGC.  The  resulting  matrix  comprises 

the  linearised  small  disturbance  equations  of  longitudinal  motion  given  in 

Ref.  (3)  Eq.  (5.1 3“ 1 y ) •  The  matrix  contains  ten  non-zero  coefficients  in 

rows  1  to  3  from  which  eleven  unknown  longitudinal  derivatives  are  to  be 

determined.  It  is  also  noted  that  two  elements  in  column  4  become  zero 

for  zero  flight  path  angle  and  that  the  denominators  in  rows  2  and  3 

contain  the  term  (2p  +  )  in  which  2p  is  almost  three  orders  of 

a 

magnitude  greater  tnan  C^  .  The  equations  are  therefore  ill-conditioned 

d 

and  in  order  to  obtain  satisfactory  estimates  for  the  eleven  unknowns,  a 

number  of  assumptions  are  proposed. 

Two  options  are  included  in  the  program.  the  first  option  can 

only  be  used  if  the  flight-path  angle  Y  is  significantly  greater  than  zero 

and  includes  the  assumption  CT^  =  ~3  Cj  for  propellor  driven  aircraft  and 

CT  =  -2  Cx  for  jet  or  rocket  powered  aircraft.  (Ref.  (p)  Section 

lV  le 

7.8).  The  remaining  ten  derivatives  are  obtained  from  the  ten  matrix 


elements.  Evaluation  of  this  method  has  shown  that  the  derivatives  C. 

L  • 

a 

and  C.  become  inaccurate  if  |Y|  <  2.0  degrees,  and  that,  the  other 

L 

q 

derivatives  are  unreliable  if  I Y |  <0.2  degrees. 


In  the  second  option,  it  is  further  assumed  that  tne  derivatives 

C  and  C,  are  related  tu  0  and  C,  directly  through  the  downwasn 

m.  L.  m  L 

u  a  q  q 

rate  dc/du  (Ref.  (3)  Section  7.10)  ana  that,  c  anu  C  are  related  to 

mm. 
q  a 

•C  and  C.  via  the  taiiplane  non-dimensional  moment  arm.  These 

Li  Li  • 

q  a 

additional  constraints  assume  that  the  contribut ions  to  tne  (  and  1 


derivatives  are  entirely  due  to  the  taij.pi.dno  lift.  Seven  ol'  tne  at. mi  wn 


derivatives  are  derived  from  the  matrix  elements  and  the  remainder  from 
the  above  assumed  relationships. 

The  second  option  provides  more  accurate  estimates  of  the 
derivatives  for  cases  where  the  assumed  constraints  are  applicable. 

5.  CONCLUSIONS 

The  six  degrees  of  freedom  dynamic  equations  of  aircraft  motion 
have  been  programmed  using  the  Advanced  Continuous  Simulation  Language  for 
use  in  aircraft  simulations.  The  air-path  axes  system  was  chosen  for  the 
integration  of  the  force  equations,  wnile  the  moment  equations  are 
integrated  in  body  axes.  Euler  angles  ana  direction  cosines  have  been 
calculated  by  use  of  quaternion  components. 

Details  of  an  example  application  have  been  provided  to 
illustrate  the  use  of  the  program  and  its  potential.  Time  histories, 
eigen  ana  jacobian  analyses  have  been  demonstrated. 
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FIG  3:  BLOCK  DIAGRAM  OF  COM8INEO  AIR-PATH  AXES/BODY  AXES 
SYSTEM  FOR  AIRCRAFT  MOTION 
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APPENDICES 

********** 

APPENDIX  1.  SDOFAP  PROCRAM  LISTING 

SDOFAP . ACSL 

PROGRAM  SDOFAP  (FLIGHT  SIMULATION  MODEL) 


"  PROGRAM  NAME  :  Six  Degrees  of  Freedom  in  Air  Path  axes.  " 
"  WRITTEN  BY  :  P.W.GIBBENS  (EOl,  A.R.L..  ABS-FW)  " 
"  COMPLETED  :  8  MAY  1985  " 
"  DESCRIPTION  THIS  PROGRAM  PRESENTS  A  SET  OF  FLIGHT  DYNAMIC  " 
"  EQUATIONS  (SIX  DEGREES  OF  FREEDOM) ,  IN  AIR-PATH  AXES,  FOR  " 
"  AIRCRAFT  SIMULATION  USING  ADVANCED  CONTINUOUS  SIMULATION  " 
"  LANGUAGE  (ACSL)  ON  THE  A.R.L.  ELXSI  6400  COMPUTER.  " 


INTEGER  ITLIM, MAXRES 

LOGICAL  BEGIN 

ARRAY  VECIN (20) . XH (3) 


INITIAL  SECTION 


INITIAL 

"»»»**»*•**  PREPARE  ALL  DATA  AND  RUNTIME  PARAMETERS  ***** 

"***  READ  IN  CONFIGURATION  AND  FLIGHT  CONDITION  DATA  ***" 
PROCEDURAL  (PHID0, ALT0, VE0K, I XX, IYY, IZZ , IXZ,MASS=) 

CALL  CAFCDI  (PHID0. ALT0, VE0K . IXX, IYY, I ZZ , IXZ , MASS) 

END  $"OF  CAFCDI  PROCEDURAL" 


"***  DEFINE  ALL  PRESET  VARIABLES  ***" 

"***  RUNTIME  PARAMETERS  ***” 

CONSTANT  $"####  RUNTIME  CONTROL  PARAMETERS  SHOULD  #### 

CONSTANT  $"####  BE  DEFINED  AT  THIS  POINT.  SEE  #### 

CONSTANT  $"####  APPLICATION  FOR  EXAMPLES.  #### 

"***  TRIMMING  ROUTINE  ITERATION  PARAMETERS 


CONSTANT 


$”####  TRIMMING  ROUTINE  PARAMETERS  SHOULD  ### 
$"####  BE  DEFINED  AT  THIS  POINT.  SEE  ### 
$”####  APPLICATION  FOR  EXAMPLES.  ««# 


■>  a  .vav  a  'A  Lva7,,a".».’^,A’.i»'.,.,.v.*.".'.’.».-  <"  ' 
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"*»*  ATMOSPHERIC  STANDARDS  AND  CONVERSION  FACTORS  ***" 
CONSTANT  RHO0  =  1.2256 


"GENERAL  CONVERSION  FACTORS,  CONSTANTS  AND  VARIABLES" 


CONSTANT 


DEGTOR= 
, KTOMPS= 
, MTONM  = 
, FTTOM  = 
,  PI  = 


0.017453 
0.514773 
0.000538 
0 . 304800 
3.141593 


RADTOD=  57.29578 
MPSTOK=  1.942602 
NMTOM  =  1858.5 
G  =  9 . 807 


...WHERE  DEGTOR=  DEGREES  TO  RADIANS 
RADTOD=  RADIANS  TO  DEGREES 
KTOMPS=  KNOTS  TO  METRES  PER  SECOND 
MPSTOK=  METRES  PER  SECOND  TO  KNOTS 
MTONM  =  METRES  TO  NAUTICAL  MILES 
NMTOM  =  NAUTICAL  MILES  TO  METRES 
FTTOM  =  FEET  TO  METRES  " 


"***  INITIAL  CONDITIONS  (FLIGHT  DATA)  **•" 


CONSTANT  DPSID0  =0.0 
,  VWN  =0.0 
,  DVE0K  =0.0 
,  PSID0  =0.0 


,  DPHID0  =0.0 
,  VWE  =0.0 
,  DX0  =0.0 
X0  =0.0 


,  DALT0  =0.0 
,  VWD  =0 . 0 
,  DY0  =0.0 
,  Y0  =0.0 


"•••  ALLOW  PARAMETER  INCREMENTING  **»" 


"####  INCREMENTS  TO  FLIGHT  CONDITION  PARAMETERS  ####" 
"####  SHOULD  BE  DEFINED  AT  THIS  POINT.  THOSE  ####" 
"###*  SHOWN  PERTAIN  TO  THE  EXAMPLE  APPLICATION.  ####" 

CONSTANT  DELI  =0.0  ,  DRPM  =0.0  ,  DXCCP=0.0  ... 

, DPL  =0.0  ,  DELS  =0.0  ,  DEL6  =0.0 


"».*  SET  CONTROL  INPUT  PARAMETERS  *«*" 


"####  REQUIRED  CONTROL  INPUT  PARAMETERS  SHOULD  ####" 
"####  BE  DEFINED  AT  THIS  POINT.  THOSE  SHOWN  ####" 
"####  PERTAIN  TO  THE  EXAMPLE  APPLICATION.  ####" 


CONSTANT  ESTART=0 . 01 , TSTART=0 . 01 
CONSTANT  EPULSE=5.0  ,TPULSE=120. 
CONSTANT  EREPET=200. ,TREPET=200. 
CONSTANT  ETAMAX=5.0  ,THRMAX=1.0 


$"**  PULSE  START  TIME  **" 
$"**  PULSE  DURATION  **" 
REPETITION  TIME  **" 
$"*•  PULSE  AMPLITUDE  **” 


PULMAX  =ETAKAX*DEGTOR 


$"*•  CONVERT  TO  RADIANS 


INCREMENT  DATA  AND  PARAMETERS 


"**•  ADD  RUNTIME  INCREMENTS  TO  FLICHT  DATA  ***" 

PSID0  =PSID0  ♦ DPS I D0 
PHID0  =PHID0  *DPHID0 
ALT0  = (ALT0  ‘DALT0) ‘FTTOM 
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VE0K  =VE|3K  +DVE0K 
X0  =XP  +DX0 
Y0  =Y0  +DY0 
GAMMAR=GAMMAR+DGAMMR 


DEFINE  INITIAL  Z  COORDINATE.  ***" 
Z0  =-ALT0 


"***  ADD  RUNTIME  INCREMENTS  TO  PARAMETERS  ***" 

VECIN (1) =ALT0  $"####  INCREMENTS  OF  FLIGHT  CONDITION  DATA  #### 
VECIN  (2)  =DRPM  $"####  ARE  ADDED  AT  SUBROUTINE  LEVEL.  #### 
VECIN (3) =DXCGP 

VECIN  (4) =DPL  $"####  INITIAL  ALTITUDE  IS  ALSO  REQUIRED  ####" 
VECIN (5) =DEL5  $"####  AT  SUBROUTINE  LEVEL.  ####" 

VECIN (6) =DEL6 

CALL  PARINC (VECIN, CGPOS.PLS) 

"***  CONVERT  FROM  EAS  TO  TAS  ***” 

PROCEDURAL  (RHO=ALT0) 

CALL  ATMOS  (ALT0.RHO) 

$"  OF  ATMOS  PROCEDURAL” 

VT0K=VE0K*SQRT  (RHO0/RHO) 

VE0  =  VE0K*KTOMPS 
VTJ3  =  VT0K*KTOMPS 


SPECIFY  SYSTEM  EXCITATION  PARAMETERS 


"***  SPECIFY  TERMINATION  CONDITION  *»*" 
CONSTANT  TSTOP  =0 . 0 


"***  SPECIFY  INDEPENDENT  VARIABLE  AS  A  PRECAUTION  ***" 

VARIABLE  TIME  =  0.0 
CINTERVAL  CINT  =  0.0 
NSTEPS  NSTP  =  0 


APPROXIMATE  TRIM  VALUES  FOR  INITIAL  CONDITIONS 


PSIR0  =  PSID0*DEGTOR 
PHIR0  =  PHID0‘DEGTOR 

PROCEDURAL  (P0 , Q0 , R0 ,  BETAR0 .  ALPHR0 ,  THETR0 .  ETAR0 
,  GAMMAR=VT0 ,  PHIR0 ,  C) 

CALL  TRAP  (VT0,PHIR0.C,GAMMAR.P0,Q0,R0.  BETAR0 
,  ALPHR0 ,  THETR0 .  ETAR0) 


$"  OF  TRAP  PROCEDURAL 


"***  CALCULATE  INERTIAL  CONSTANTS  **** 

Cf)  =  ((IXX*IZZ)-(IXZ*IXZ)) 

Cl  =  IZZ/C0 
C2  =  IXZ/C0 
C3  -  C2*(IXX-IYY*IZZ) 

C4  =  ( (IYY-IZZ) *C1) - (IXZ*C2) 

C5  =  1.0/IYY 
C6  =  C5*IXZ 
C7  =  C5* (IZZ-IXX) 

C8  =  IXX/C0 

C9  =  ( (IXX-IYY) *C8) ♦ (IXZ*C2) 


INITIALIZE  QUATERNIONS 


THETD0=  THETR0*RADTOD 


1 p  . .  CONTINUE 

"**NOTE: -THE  EULER  ANGLE (S)  ARE  HALVED  WHEN  CALCULATING... 
THE  QUATERNIONS" 

CALL  QUATNS (PSIR0, THETR0, PHIR0, TAU00 , TAU10, TAU20 , TAU30) 


END  $"OF  INITIAL" 


DYNAMIC  SECTION 


DYNAMIC 

”***  THE  TRANSITION  FROM  INITIAL  TO  DYNAMIC  TRANSFERS  ***” 
"»**  ALL  INITIAL  CONDITIONS  TO  THE  STATE  VARIABLES  AND  ***" 
"»**  EVALUATES  THE  CODE  IN  THE  DERIVATIVE  SECTION  ONCE.  ***” 


IF (BEGIN)  GO  TO  20 


TRIM  AIRCRAFT  WITH  USER  SUPPLIED  SUBROUTINE 


"####  EXAMPLE  OF  TRIMMING  IN  RECTILINEAR  ####" 

"####  FLICHT  USING  SUBROUTINE  POWIT.  ####" 

XH (1)  =THETR0  $"####  OTHER  TRIM  CONDITIONS  MAY  BE  #### 

XH (2)  =ALPHR0  $”####  OBTAINED  BY  DEFINING  APPROPRIATE  4M *## 

XH (3)  =ETAR0  $"####  TRIM  STATES  AND  ASSOCIATED  TRIM  #### 

$"####  CONTROLS.  EC.  LEVEL  BANKED  TURN.  #### 


PROCEDURAL (XH=XH, MAXRES . ITLIM. ERRMAX) 

CALL  POWITfXH, MAXRES. ERRMAX. ITLIM) 

THETR0  -XH ( 1 ) 

ALPHR0  =XH ( 2 ) 

ETAR0  =XH (3) 


END 


$  "OF  POWIT  PROCEDURAL 


BEGIN=.TRUE. 
GO  TO  10 

20  . .  CONTINUE 


DERIVATIVE  SECTION 


DERIVATIVE 


CALCULATE  CONTROL  INPUTS 


"####  EXAMPLE  OF  CONTROL  INPUTS  USINC  THE  ACSL  PULSE  ####" 
"####  GENERATOR,  OTHER  INPUT  FORMS  ARE  AVAILABLE.  ####" 

PROCEDURAL (ETAR . PLS= 

ETAR0, DELPL, ESTART, EREPET, EPULSE , TSTART, TREPET, TPULSE) 

DELETA=PULMAX*PULSE (ESTART, EREPET, EPULSE) 

DELPL  =THRMAX*PULSE (TSTART, TREPET, TPULSE) 

ETAR=ETAR0*DELETA 

CALL  CONTROLS (ETAR . DELPL . PLS) 

ETAD=ETAR*RADTOD 

END  $  "OF  CONTROLS  PROCEDURAL" 


"***•»***«  CALCULATE  QUATERNIONS;  NORMALIZE,  AND  THEN  *»** 
. .  DETERMINE  THE  DIRECTION  COSINES  **•* 

PROCEDURAL  (LI, L2.L3, Ml. M2, M3, N1,N2 , N3=TAU0, TAU1, TAU2, TAU3) 


NORMALIZE  QUATERNIONS  ***" 

TAUN  =  SQRT ( (TAU0*TAU0) ♦ (TAU1 *TAU1) ♦ (TAU2*TAU2) ♦ (TAU3‘TAU3) ) 

TAU0N  =  TAU0/TAUN 

TAU1N  =  TAU1/TAUN 

TAU2N  =  TAU2/TAUN 

TAU3N  =  TAU3/TAUN 

"•**  CALCULATE  THE  DIRECTIONAL  COSINES  (LI  ->  N3)  ***" 

LI  =  ( ( (TAU0N*TAU0N) ♦ (TAU1N*TAU1N) ) *2.0) -1.0 
L2  =  ( (TAU1N*TAU2N) ♦ (TAU0N*TAU3N) ) *2 .0 

L3  =  ( (TAU1N*TAU3N) - (TAU0N‘TAU2N) )  *2 .0 

Ml  =  ( (TAU1N* TAU2N) - (TAU0N*TAU3N) ) *2 .0 

M2  =  (  ( (TAU0N»TAU0N) ♦ (TAU2N*TAU2N)  )  *2 .0) -1 .0 
M3  =  ( (TAU2N*TAU3N) ♦ (TAU0N*TAU1N) ) *2 .0 

N1  =  (  (TAU1N*TAU3N) * (TAU0N*  TAU2N) ) *  2 . 0 

N2  =  ( (TAU2N*TAU3N) - (TAU0N * TAU1N)  )  *2 .0 

N3  =  ( ( (TAU0N*TAU0N) ♦ (TAU3N*TAU3N) ) *2 .0) -1 .0 


END  $"OF  QUATERNION  PROCEDURAL 


"***  CALCULATE  ATMOSPHERIC  CONDITIONS  FOR  CURRENT  ALTITUDE  *** 


CALL  ATMOS  (-Z,RHO)  $"####  ATMOSPHERICS  SUBROUTINE  MUST  ####" 

"####  BE  SUPPLIED  BY  THE  USER.  ####" 


CALCULATE  LONGITUDINAL  AND  LATERAL  AERODYNAMIC  FORCES  ***" 
"“****“  IN  STABILITY  AXES  AND  MOMENTS  IN  BODY  AXES  *««*****»" 


PROCEDURAL  (XA, YA, ZA, LA, MA, NA=VT , ALPHAR , BETAR , P , Q , R , PHIR , THETAR . DALPHR) 

"####  THE  AERO  SUBROUTINE  SUPPLIES  THE  AERODYNAMIC  FORCES  ####" 
"####  IN  STABILITY  AXES.  THIS  MUST  BE  SUPPLIED  BY  USER.  ####" 

CALL  AERO  (XA, YA, ZA, LA, MA, NA, VT, ALPHAR, BETAR, P, Q, R, 

PHIR , THETAR , DALPHR) 

END  $"OF  AERO  PROCEDURAL" 

PROCEDURAL  (XP, YP, ZP, LP , MP, NP=VT, ALPHAR, BETAR, P,Q,R, PHIR, THETAR) 

"####  THE  PROP  SUBROUTINE  SUPPLIES  THE  PROPULSIVE  FORCES  ####" 
"####  IN  BODY  AXES.  THIS  MUST  BE  SUPPLIED  BY  USER.  ####" 

CALL  PROP  (XP,YP,ZP.LP,MP,NP,VT, ALPHAR, BETAR. P,Q,R,  ... 

PHIR, THETAR) 

END  $"OF  PROP  PROCEDURAL" 

SA  =  SIN (ALPHAR) 

CA  =  COS (ALPHAR) 

SB  =  SIN  (BETAR) 

CB  =  COS (BETAR) 


"***  CALCULATE  GRAVITY  FORCE  COMPONENTS  ->  STABILITY  AXES  ***" 

FGX  =  L3*MASS*G*CA  ♦  N3*MASS*G*SA 
FGY  =  M3*MASS*C 

FGZ  =-L3*MASS*G*SA  ♦  N3*MASS*G*CA 


”***  CALCULATE  THE  TOTAL  FORCES  (FXS , FYS, FZS)  ***" 
"***  ->  STABILITY  AXES  . 

FXS  =  XA  ♦  XP*CA  +  ZP*SA  ♦  FGX 
FYS  =  YA  *  YP  +  FCY 
FZS  =  ZA  -  XP‘SA  ♦  ZP*CA  ♦  FGZ 


"*»*  CALCULATE  THE  TOTAL  MOMENTS  (L,M,N)  ->  BODY  AXES  ***" 

L  =  (LA‘CA)  -  (NA*SA)  ♦  LP 
M  =  MA  ♦  MP 
N  =  (LA*SA)  ♦  (NA‘CA)  ♦  NP 


»  #  *  * 
"tit 


CALCULATE  THE  TOTAL  FORCES  (FXW, FYW, FZW) 
->  AIR-PATH  AXES 


***  CALCULATE  FORCES  IN  BODY  AXES  ***" 


FXB  =  XA*CA  -  ZA*SA  ♦  XP 

FYB  =  YA  ♦  YP 

FZB  =  XA*SA  ♦  ZA*CA  ♦  ZP 


"***  CALCULATE  LINEAR  ACCELERATIONS  ***" 

AX  =  FXB/ (MASS *G) 

AY  =  FYB/(MASS«G) 

AZ  =  FZB/(MASS*G) 


"***  CALCULATE  NORMAL  ACCELERATION  ***" 
AN  =  -AZ 


"***  RESOLVE  BODY  AXIS  ANGULAR  RATES  INTO  STABILITY  AXES  ***" 

PS  =  P*CA  +  R*SA 
RS  =  -P*SA  ♦  R*CA 


CALCULATE  LINEAR  DERIVATIVES  ->  AIR  PATH  AXES  ***“*" 

"***  DECLARE  DALPHR  AS  AN  IMPLICIT  VARIABLE  ***" 

DALPHR  =  IMPL  (Q, 0 . 001 , 30, EF , 

(Q*C8-PS«SB*FZW/MASS/VT) /CB , 0  -  0001) 

PROCEDURAL  (=EF) 

CALL  IMP  (EF) 

END  $"OF  DALPHR  PROCEDURAL" 


DBETAR  =  FYW/MASS/VT-RS 
DVT  =  FXW/MASS 


"***  CALCULATE  ANGULAR  ACCELERATIONS  ->  BODY  AXES  ***” 

PDOT  =  (L*C1) ♦ (N*C2) ♦ ( ( (P*C3) ♦ (R*C4) ) *Q) 

QDOT  =  (M*C5) ♦ ( ( (R*R) - (P*P) ) *C6) ♦ (R*P*C7) 

RDOT  =  (N*C8)  ♦  (L‘C2)  ♦  (  ( (P*C9)  -  (R*C3) )  *Q) 


.  m 
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DETERMINE  QUATERNION  RATES  AND  EULER  ANGLES.  “******" 


"***  CALCULATE  THE  QUATERNION  RATES  *»*" 

TAU0DT=- ( ( (TAU1N‘P) ♦ (TAU2N*Q) ♦ (TAU3N«R) ) *0 . 5) 
TAU1DT=  ( (TAU0N*P) - (TAU3N‘Q) ♦ (TAU2N‘R) ) *0 . 5 
TAU2DT=  ( (TAU3N*P) ♦ (TAU0N*Q) - (TAU1N*R) ) *0. 5 
TAU3DT=- ( ( (TAU2N*P) - (TAUIN'Q) - (TAU0N*R) ) *0. 5) 


"»**  CALCULATE  THE  NEW  EULER  ANGLES 


PROCEDURAL (PSI D, THETAD, PHI D=L1, L2. L3. M3, N3) 


— 


/vvvvyy.*:  l- v  >  v 


"***  CALCULATE  HEADING  (OR  YAW)  ANGLE  (PSI) 

"***  IN  RANGE  (9 (NORTH)  TO  360  DEG.  ***" 

PSIR  =  ATAN (L2/L1) 

IF  (PSIR. LT. 0.0)  PSIR  =  PSIR*PI*2 
PSID  =  PSIR*RADTOD 


"*»*  CALCULATE  ANGLE  OF  PITCH  (THETA)  ***" 
"***  IN  RANGE  */-9p  DEGREES  ***" 

THETAR=  ATAN (-L3/SQRT ( (L1*L1) ♦ (L2*L2) ) ) 
THETAD=  THETAR‘RADTOD 


"**•  CALCULATE  BANK  (OR  ROLL)  ANGLE  (PHI)  IN  RANGE  ***" 
"***  ♦/-18{J  DEGREES , WHERE  0  DEG.  INDICATES  WINGS  LEVEL  ***" 

PHIR  =  ATAN (M3/N3) 

PHID  =  PHIR* RADTOD 

END  9 "OF  EULER  PROCEDURAL" 


TRAJECTORY 


"***  RESOLVE  FLIGHT  PATH  AXES  VELOCITIES  INTO  BODY  ***" 

"***  COMPONENTS  FOR  CALCULATION  OF  DISTANCES  IN  EARTH  AXES  ***" 

UB  =  VT*CA*CB 
VB  =  VT*SB 
WB  =  VT*SA*CB 


"*•*  CALCULATE  EARTH  AXES  VELOCITIES  FROM  BODY  AXES  ***" 
"***  VELOCITIES  USING  DIRECTION  COSINES  AND  WIND  VELOCITIES  ***" 

VNE  =  L1*UB+M1*VB+N1*WB-VWN 
VEE  =  L2*UB*M2*VB*N2*WB-VWE 
VDE  =  L3*UB+M3*VB+N3*WB- VWD 


PROCEDURAL (VGRKT , GAMMAD , CHID=VNE , VE  E . VDE ) 


"»**  CALCULATE  THE  VELOCITY  OVER  THE  GROUND,  *  *  *  " 
"»»*  VGR  (IE.  GROUND  SPEED)  ***" 

VCR  =  SQRT((VNE*VNE)+(VEE*VEE)) 

VCRKT  =  VGR*MPSTOK 


'***  CALCULATE  FLICHT-PATH  ANGLE,  GAMMA  IN  RANGE  */-9?  DEG  *•*" 

GAMMAR  =  THETAR-ALPHAR 
GAMMAD  =  CAMMAR ‘ RADTOD 


'**•  CALCULATE  ANGLE  OF  CLIMB,  LAMBDA  IN  RANGE  DEG  *»*" 

LAMDAR  =  ATAN (- VDE /VCR) 

LAMDA D  =  LAMDAR* RADTOD 
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CALCULATE  ANGLE  OF  TRACK,  CHI,  ***" 
"***  IN  RANGE  0 (NORTH)  TO  360  DEG.  ***” 

CHIR  =  ATAN (VEE/VNE) 

CHID  =  CHIR*RADTOD 

END  $"OF  TRAJECTORY  PROCEDURAL" 


"***  CALCULATE  THE  RANGE  (NOTE:  IN  METRES)  ***" 
RANGE  =  SQRT ( (X*X) + (Y*Y) ) 


INTEGRATION  OF  SYSTEM  STATE  EQUATIONS 


ALPHAR 

INTVC (DALPHR , ALPKR0) 

BE  TAR 

= 

INTVC (DBETAR, BETAR0) 

VT 

— 

INTVC (DVT, VT0) 

P 

INTVC (PDOT.P0) 

$"ROLL  RATE" 

Q 

= 

INTVC (QDOT.Q0) 

$"PITCH  RATE" 

R 

s 

INTVC (RDOT.R0) 

$"YAW  RATE" 

TAU0 

= 

INTVC (TAU0DT. TAU00) 

$"QUATERNION  TERMS" 

TAU1 

= 

I NTVC ( TAU IDT , TAU 10) 

TAU2 

= 

INTVC (TAU2DT, TAU20) 

TAU3 

I NTVC (TAU3DT , TAU  30) 

Z 

= 

I NTEG ( VDE , Z0) 

$"Z-POSITIONAL  CO-ORDINATE" 

X 

= 

INTEG(VEE,X0) 

$"X- POSITIONAL  CO-ORDINATE" 

Y 

= 

INTEC (VNE , Y0) 

$"Y-POSITIONAL  CO-ORDINATE" 

END  $"OF  DERIVATIVE" 

VE 

= 

VT*SQRT (RHO/RHO0) 

$ "EQUIVALENT  AIRSPEED 

(M/S) " 

VEK 

= 

VE*MPSTOK 

$"EQUIVALENT  AIRSPEED 

(KTS) " 

ALPHAD 

= 

ALPHAR  *  RADTOD 

5 "ANGLE  OF  ATTACK 

(DEG) " 

ALT 

=  ■ 

-Z/FTTOM 

$ "ALTITUDE 

(FT.) " 

EXPRESS  TERMINATION  CONDITION 


TERMT (TIME . GE . TSTOP) 


END  $"OF  DYNAMIC" 


NOTE:  THIS  LISTING  SHOULD  BE  USED  IN  CONJUNCTION  WITH  ***" 
"***  THE  ADVANCED  CONTINUOUS  SIMULATION  LANGUAGE  (ACSL)  USER  ***" 
"*»*  GUIDE/REFERENCE  MANUAL. 


END  $"OF  PROGRAM 


n  n  •»  no 


>  V  >  V 
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SUBROUTINE  QUATNS (PSIR , THETAR , PHIR , TAU0, TAU1 , TAU2 . TAU3) 

C  "***  CALCULATES  NORMALIZED  QUATERNIONS  »**" 

SPSI  =  SIN  (PSIR*0 . 5) 

CPSI  =  COS(PSIR*0.5) 

STKETA  =  SIN (THETAR *0.5) 

CTHETA  =  COS (THETAR*0 . 5) 

SPHI  =  SIN  (PHIR*0 . 5) 

CPHI  =  COS (PHIR*0 . 5) 


•  V-  -o 
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C  "***  CALCULATE  THE  QUATERNIONS  FROM  THE  EULER  ANGLE  TERMS  ***" 


TAU0  =  (CPHI*CTHETA*CPSI) ♦ (SPHI *STHETA*SPSI) 
TAU1  =  (SPHI *CTHETA*CPSI )- (CPHI* STHETA* SPSI) 
TAU2  =  (CPHI* STHETA* CPS I) ♦ (SPHI *CTHETA*SPSI) 
TAU3  =  (CPHI *CTHETA*SPSI) - (SPHI*STHETA*CPSI) 


C  "***  NORMALIZE  INITIAL  QUATERNIONS  ***" 

TAUN  =  SQRT ( (TAU0*TAU0) ♦ (TAU1*TAU1) ♦ (TAU2*TAU2) 
1  ♦ (TAU3*TAU3) ) 

TAU0  =  TAU0/TAUN 
TAU1  =  TAU1/TAUN 
TAU2  =  TAU2/TAUN 
TAU3  =  TAU3/TAUN 


RETURN 


END 


SUBROUTINE  EVAL  (XXXX. XXXRES , MAXR) 

"***  THIS  SUBROUTINE  IS  CALLED  FROM  THE  TRIMMING  SUBROUTINE  ***" 
"*»*  TO  GAIN  ACCESS  TO  THE  DERIVATIVE  SECTION.  ***" 

DIMENSION  XXXX (3) , XXXRES (3) 

’’***  THE  "9"  CONTROL  CHARACTER  MAKES  MAIN  PROGRAM  ***" 

"***  VARIABLES  AVAILABLE  TO  THIS  SUBROUTINE.  ***" 

INTEGER  I 
1  =  1 


THETAR  =  (XXXX(l)) 

ALPHAR  =  (XXXX (2)) 

ETAR0  =  (XXXX (3)) 

CALL  QUATNS (PSIR0 , THETAR . PHIR0, TAU0, TAU1 , TAU2 . TAU3) 
CALL  ZZDERV(I) 


xxxres  (i)  =  (dvt) 

XXXRES  (2)  =  (DALPHR)  >Vt- 

XXXRES  (3)  =  (QDOT)  'vS'v  -V 


RETURN 


SUBROUTINE  IMP(EF) 


NOTIFIES  FAILURE  OF  IMPLICIT  ITERATION  ROUTINE 


IF  (EF  .EQ.  0.0)  GOTO  89 
WRITE  (6,99) 

FORMAT (//27HIMPLICIT  PROCEDURAL  FAILED  //) 


RETURN 
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FTSUBS.F 


SUBROUTINE  CAFCDI (PHIN.HNM, VN, IXX, I YY, IZZ , 1X2, MASS) 

IMPLICIT  REAL (A-Z) 

INCLUDE  'FTPAR.F' 

INTEGER  I,J 


C  ******  ROUTINE  FOR  THE  INPUT  OF  CONFIGURATION  AND  ****** 

C  ***************  FLIGHT  CONDITION  DATA.  ***************** 


«V  *fct  •* 
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OPEN  (UNIT=1 , FILE= 'FTSUD2 .IN', STATUS= ' OLD 1 ) 
OPEN  (UNIT=7, FILE= 'FTGRAF . VPP ' , STATUS= ’OLD ' ) 


C  ***  READ  IN  CONFIGURATION  AND  FLIGHT  CONDITION  DATA  *** 

C  ####  THE  VARIABLES  SHOWN  PERTAIN  TO  THE  EXAMPLE  APPLICATION  #### 

C  ####  DATA  READ  IN  IS  AVAILABLE  AT  SUBROUTINE  LEVEL  ONLY.  «««# 


READ ( 1 ,*) TTOT , TSAMP , NH , HN , DELH , NV , VN , DELV . NPHI . PHI N 

READ ( 1 ,  * ) DELPHI , NPL , PL?, DELPL . NRPM , RPM , DELRPM , WECHT . NXCG . XCGP 

READ (1, * ) DELXCG, ZCG, (PEFF(I) ,1=1,7) .XT 

READ ( 1 , * ) ZT , XTH , ZTH , THSET , MAXEP , NBLAD , PDI A . WSET . TPSET . ETAG 
READ ( 1 , * ) ACMASS , ACIXX , ACI YY, ACI ZZ , ACIXZ , SW, CW. BW, ST, CT 
READ ( 1 ,*) CETA , XP , ZP , CWP , BFW , BTAI L , CL? , CLAL . CD? . CDAL 
READ  ( 1 ,  * )  CM? .  CMAL ,  EPS? ,  EPSAL ,  QTOQ . A? , A1 . A2 . A3 . B? 

READ ( 1 , * ) B1 , B2 , B3 , CD?T , CDLT , CMT? , CMQW , CLP , CLXI ,  KWI NG 

READ ( 1 , * ) KFUSE . NPSF I . NPSFR , MPSF I , MPSFR , CDB , CTR . BL? . WTR , TAPERF 

READ  ( 1 .  * )  XWC ,  ZWC ,  S?S ,  XQARTC 

READ  (7, »)  (FPROP(I.l) ,FPROP(I,2) , 1  =  1,2?) 

READ  (7, *)  ((TOP (I, J) , J=l, 33) . 1  =  1, 2?) 

READ (7,*) ((ETP(I.J) , J=1 , 6?) . 1=1 , 21) 

111  READ (7.*) (CYPROP (1,1) ,CYPR0P(I,2) ,CYPROP(I,3) ,CYPR0P(I,4) ,1=1,2?) 
READ (7, »)  ( (DELEPS (I , J) . J=l. 1?) ,1=1,8) 

C  ####  DATA  WHICH  IS  REQUIRED  IN  THE  MAIN  PROGRAM  MUST  #### 

C  ####  RENAMED  BEFORE  BEING  PASSED  AS  ARGUMENTS.  «««# 

MASS=ACMASS 
IXX=ACIXX 
IYY=ACIYY 
I Z7=ACI ZZ 
IXZ=ACIXZ 

HNM=HN* . 3?48 

CLOSE  (UNIT=1) 

CLOSE (UNIT=7) 


RETURN 
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SUBROUTINE  PARINC  (VECIN, CGPOS.PLS) 

"*“  SUBROUTINE  TO  ALLOW  RUNTIME  PARAMETER  MODIFICATIONS  ***" 

IMPLICIT  REAL (A-Z) 

DIMENSION  VECIN (6) 

INCLUDE  ' FTPAR . F ' 

####  THE  VARIABLES  SHOWN  PERTAIN  TO  THE  EXAMPLE  APPLICATION.  #### 
####  THEY  DEMONSTRATE  THE  METHOD  OF  PARAMETER  INCREMENTING.  #### 

ALT0=  VECIN (1) 

RPM  =RPM  + VECIN (2) 

XCGP=XCGP+VECIN (3) 

PL0  =PL0  ♦ VECIN  (4) 

PAR5=PAR5+VECIN (5) 

PAR6=PAR6+ VECIN (6) 

XCG=XCGP  *  CW/ 100 . +XQARTC-CW/4 . 

####  DATA  WHICH  IS  REQUIRED  IN  THE  MAIN  PROGRAM  MUST  #### 

####  RENAMED  BEFORE  BEING  PASSED  AS  ARGUMENTS.  #### 

CGPOS=XCGP 

PLS=PL0 
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RETURN 

END 

SUBROUTINE  TRAP (VT0 , PHIR0 . G . GAMMAR . P0 , Q0 . R0 , BE TAR 0 


1  , ALPHR0, THETR0, ETAR0) 

C  ********  GIVES  INITIAL  APPROXIMATION  FOR  TRIM  CONDITION 


C  ####  THE  APPROXIMATION  SHOWN  PERTAINS  TO  THE  #### 

C  ####  EXAMPLE  APPLICATION.  THE  APPROXIMATION  TO  #### 

C  ##*#  BE  USED  IS  AT  THE  USERS  DISCRETION.  #### 

IMPLICIT  REAL (A-Z) 

INCLUDE  ' FTPAR. F ' 

PI  =  3.141593 
ETA00=  2.2 

KK  =41 .0* (XCGP/100 .-0.4) 

QD  =(.5*1 .2256*VT0*VT0) 

P0  =0. 

Q0  =G/VT0*TAN (PHIR0) *SIN (PHIR0) 

R0  =G/VT0*SIN (PHIR0) 

AN  =  1/COS (PHIR0) 

BETAR0  =  0. 


ALPHAW  =  ACMAS  S  *  G/ (CLAL  *  QD  *  SW) 
ALPHR0  =  ALPHAW  -  WSET 
CLWB  =  CL0  ♦  (CLAL ‘ALPHAW) 
CDWB  =  CD0  ♦  (CDAL*CLWB*  *  2) 

DRAG  =  CDWB  *  QD  *  SW 

CAMMAR  =  ASIN ( -DRAG/ (ACMASS*G) ) 
THETR0  =  GAMMAR  ♦  ALPHR0 
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ETAp  =  (ETA00*KK*CLWB) *PI/180. 
ETAR0  =  ETAp 


RETURN 

END 


SUBROUTINE  ATMOS  (H,RO) 

IMPLICIT  REAL (A-Z) 

COMMON/ARGS2/HITE , DENS 

C  ****  THIS  ROUTINE  GIVES  THE  DENSITY  AT  A  **** 

C  ****  GIVEN  ALTITUDE  FOR  ISA  CONDITIONS.  **** 

C  ####  MORE  COMPLICATED  ATMOSPHERE  ROUTINES  MAY  BE  USED.  #### 

C  ####  THE  SUBROUTINE  TO  BE  USED  IS  AT  THE  USERS  DISCRETION.  #### 


HI TE=H 

DENS=1 . 2256*  (1-0.02256*HITE/1000.0) **4.2561 
RO  =DENS 


RETURN 

END 


SUBROUTINE  AERO  (XA. YA, ZA, LA, MA, NA, VT , ALPHAR , BETAR 


1  . PP , QQ , RR . PHIR , THE TAR) 

C  ********  CALCULATES  TOTAL  AERODYNAMIC  FORCES  AND  MOMENTS 


IMPLICIT  REAL (A-Z) 
INCLUDE  ’FTPAR.F' 
DIMENSION  X (20) 


C  ####  THE  VECTOR  X  IS  USED  TO  PASS  STATE  VARIABLES  TO  #### 

C  ####  ANY  SUBROUTINES  CALLED  FROM  AERO.  SEE  THE  EXAMPLE  #### 

C  ####  APPLICATION  FOR  DETAILS  OF  ITS  USE.  #### 

X(V)  =  VT 
X (ALPHA)  =  ALPHAR 
X(Q)  =  QQ 
X(P)  =  PP 
X  (H)  =  HITE 
X  (PHI )  =  PHIR 
X (THETA)  =  THETA R 

CALL  ATMOS (HI TE . RHO) 

QD= . 5*RHO*VT* • 2 

C  ********  DETERMINE  INDIVIDUAL  AERODYNAMIC  CONTRIBUTIONS.  ** 


C  ####  ANY  SUBROUTINE  CALLS  OR  OTHER  CALCULATIONS  «««« 
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*«««  TO  DETERMINE  INDIVIDUAL  AERODYNAMIC  FORCE  *##* 
*#**  CONTRIBUTIONS  SHOULD  BE  ENTERED  HERE.  *#** 
«#«#  SEE  APPLICATION  FOR  EXAMPLES.  *##* 


C 


*******  CALCULATE  TOTAL  AERODYNAMIC  FORCES  AND  MOMENTS. 

*#**  TOTAL  FORCES  AND  MOMENTS  SHOULD  BE  DETERMINED  **## 
####  HERE  BY  ADDITION  OF  INDIVIDUAL  CONTRIBUTIONS.  #### 
«##*  SEE  APPLICATION  FOR  EXAMPLES.  ###« 

XA  = 

YA  = 

ZA  = 


LA  = 
MA  = 
NA  = 


RETURN 

END 


SUBROUTI NE  PROP (XP . YP , ZP . LP , MP , NP , VT , ALPHAR , BETAR , 

1  PP.QQ,RR,PHIR,THETAR) 

********  CALCULATES  TOTAL  PROPULSIVE  FORCES  AND  MOMENTS 

IMPLICIT  REAL (A- Z) 

INCLUDE  'FTPAR.F' 

DIMENSION  X(20) 


####  THE  VECTOR  X  IS  USED  TO  PASS  STATE  VARIABLES  TO  #««# 
####  ANY  SUBROUTINES  CALLED  FROM  PROP.  SEE  THE  EXAMPLE  «#«* 
####  APPLICATION  FOR  DETAILS  OF  ITS  USE.  #### 

X(V)  =  VT 
X (ALPHA)  =  ALPHAR 
X(Q)  =  QQ 
X(P)  =  PP 
X (H)  =  HITE 
X (PHI)  =  PHIR 
X  (THETA)  =  THE  TAR 

CALL  ATMOS (HI TE.RHO) 

QD= . 5*RHO*VT*  *2 


DETERMINE  INDIVIDUAL  PROPULSIVE  CONTRIBUTIONS. 


####  ANY  SUBROUTINE  CALLS  OR  OTHER  CALCULATIONS  #### 
####  TO  DETERMINE  PROPULSIVE  FORCES  SHOULD  BE  #### 
####  ENTERED  HERE.  #### 
####  SEE  APPLICATION  FOR  EXAMPLES.  #### 


*******  CALCULATE  TOTAL  PROPULSIVE  FORCES  AND  MOMENTS. 

####  TOTAL  FORCES  AND  MOMENTS  SHOULD  BE  DETERMINED  #### 
«*«*  HERE  BY  ADDITION  OF  INDIVIDUAL  CONTRIBUTIONS.  #### 
####  SEE  APPLICATION  FOR  EXAMPLES.  #### 
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RETURN 

END 


SUBROUTINE  CONTROLS (ETAR, DELPL , PLS) 

INCLUDE  ' FTPAR.F' 

***  THIS  SUBROUTINE  IS  USED  TO  INCREMENT  OR  OTHERWISE  *** 

*♦*  MODIFY  CONTROL  INPUTS  AS  REQUIRED.  *** 

####  THE  CODINC  SHOWN  RELATES  TO  THE  EXAMPLE  APPLICATION.  #### 

ETA  =  ETAR 
PL  =  PL0*DELPL 
PLS  =  PL 

RETURN 
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SDOFAP . ACSL 


PROCRAM  SDOFAP  (FLICHT  SIMULATION  MODEL) 


"  PROGRAM  NAME  :  Six  Degrees  of  Freedom  in  Air  Path  axes. 

"  WRITTEN  BY  :  P.W.  CIBBENS  (EOl,  ABS-FW,  ARL) 

"  COMPLETED  8  MAY  1985 

H 

"  MODIFIED  BY  R.H.  PERRIN  (EOl.  ABS-RW,  ARL)  TO  ALLOW  FOR 
"  A  VERY  SMALL  OR  ZERO  CLIMB  ANGLE  -  OCTOBER  1985 

n 

"  DESCRIPTION  :  THIS  PROGRAM  PRESENTS  A  SET  OF  FLIGHT  DYNAMIC 
”  EQUATIONS  (SIX  DEGREES  OF  FREEDOM).  IN  AIR-PATH  AXES,  FOR 
"  AIRCRAFT  SIMULATION  USING  ADVANCED  CONTINUOUS  SIMULATION 
"  LANGUAGE  (ACSL)  ON  THE  ARL  ELXSI  6400  COMPUTER. 


INTEGER  ITLIM, MAXRES 

LOGICAL  BEGIN 

ARRAY  VECIN  (20) , XH (3) 


INITIAL  SECTION 


INITIAL 

«»**»«*****  PREPARE  ALL  DATA  AND  RUNTIME  PARAMETERS  ***** 

"*»*  READ  IN  CONFIGURATION  AND  FLICHT  CONDITION  DATA  ***" 
PROCEDURAL  (PHID0 , ALT0 ,VE0K.IXX,IYY,IZZ,IXZ, MASS=) 

CALL  CAFCDI  (PHID0, ALT0, VE0K, I XX, IYY, IZZ, IXZ.MASS) 

END  $"OF  CAFCDI  PROCEDURAL" 


"***  DEFINE  ALL  PRESET  VARIABLES  ***" 
"»**  DEFINE  RUNTIME  PARAMETERS  ***" 


CONSTANT  TSTP=5 
CONSTANT  HI  1=0. 
CONSTANT  LO1=0 . 


.  PLS=0 .  , CGPOS=25 . 1 

,  HI  2=0 .  , HI  3=0 . 

,  LO2=0 .  , LO3=0 . 


***  POWIT  TRIMMING  ROUTINE  ITERATION  PARAMETERS  •»•" 


CONSTANT  MAXRES  =3 

.ERRMAX  =0.000000 1 


.ITLIM  =100 


"***  ATMOSPHERIC  STANDARDS  AND  CONVERSION  FACTORS  ***" 
CONSTANT  RHO0  =  1.2256 


"***  GENERAL  CONVERSION  FACTORS,  CONSTANTS  AND  VARIABLES 


CONSTANT 


DEGTOR= 
, KTOMPS= 
, MTONM  = 
, FTTOM  = 
,  PI  = 


0.017453 
0.514773 
0.000538 
0 . 304800 
3.141593 


RADTOD=  57.29578 
MPSTOK=  1.942602 
NMTOM  =  1858.5 
G  =  9 . 807 


". . .WHERE 


DEGTOR=  DEGREES  TO  RADIANS 
RADTOD=  RADIANS  TO  DECREES 
KTOMPS=  KNOTS  TO  METRES  PER  SECOND 
MPSTOK=  METRES  PER  SECOND  TO  KNOTS 
MTONM  =  METRES  TO  NAUTICAL  MILES 
NMTOM  =  NAUTICAL  MILES  TO  METRES 
FTTOM  =  FEET  TO  METRES  " 


"***  INITIAL  CONDITIONS  (FLIGHT  DATA)  **»" 


CONSTANT 

DPSID0  =0.0 

, DPHID0  =0.0 

, DALT0  =0.0 

,  VWN  =0 . 0 

,  VWE  =0.0 

,  VWD  =0.0 

.  DVE0K  =0.0 

,  DX0  =0.0 

,  DY0  =0.0 

,  PSID0  =0.0 

,  X0  =0.0 

,  Y0  =0.0 

"***  ALLOW  PARAMETER  INCREMENTING  ***" 

CONSTANT 

DELI  =0.0 

,  DRPM  =0.0 

,  DXCGP=0 . 0  ... 

. DPL  =0.0 

.  DEL5  =0.0 

,  DEL6  =0.0 

”***  SET 

CONTROL  INPUT 

PARAMETERS  ***" 

CONSTANT  ESTART=0.01,TSTART=0.01  $"**  PULSE  START  TIME  ** 
CONSTANT  EPULSE=5.0  ,TPULSE=120.  SH‘*  PULSE  DURATION  ** 
CONSTANT  EREPET=200. ,TREPET=200.  $"**  REPETITION  TIME  ** 
CONSTANT  ETAMAX=5 . 0  , THRMAX=1 . 0  $"**  PULSE  AMPLITUDE  ** 

PULMAX  =E TAMAX * DE GTOR  $"***  CONVERT  TO  RADIANS 


"***  ADD  RUNTIME  INCREMENTS  TO  FLIGHT  DATA  ***" 

PSID0  =PSID0  +DPSID0 
PHID0  =PHID0  +DPHID0 
ALT0  = (ALT0  *DALT0) *FTTOM 
VE0K  =VE0K  *DVE0K 
X0  =X0  +DX0 
Y0  =Y0  +DY0 
GAMMAR  =GAMMAR *  DCAMMR 


"  *  *  * 


DEFINE  INITIAL  Z  COORDINATE. 


W.V^'j 


I 


"•**  ADD  RUNTIME  INCREMENTS  TO  PARAMETERS  “*" 

VECIN  (1) =ALT0 
VECIN (2) =DRPM 
VECIN (3) =DXCGP 
VECIN  (4) =DPL 
VECIN  (5) =DEL5 
VECIN (6) =DEL6 

CALL  PARINC  (VECIN, CGPOS,PLS) 

"***  CONVERT  FROM  EAS  TO  TAS  ***” 

PROCEDURAL  (RHO=ALT0) 

CALL  ATMOS  (ALT0,RHO) 

END  $"  OF  ATMOS  PROCEDURAL” 

VT0K=VE0K*SQRT  (RHO0/RHO) 

VE0  =  VE0K*KTOMPS 
VT0  =  VT0K  *  KTOMPS 

. .  SPECIFY  SYSTEM  EXCITATION  PARAMETERS 

"*«*  SPECIFY  TERMINATION  CONDITION  ***" 

CONSTANT  TSTOP  =0.0 


"***  SPECIFY  INDEPENDENT  VARIABLE  AS  A  PRECAUTION  ***" 

VARIABLE  TIME  =0.0 
CINTERVAL  CINT  =0.0 
NSTEPS  NSTP  =  0 


"***  SET  INITIAL  APPROXIMATION  FOR  INITIAL  CONDITIONS  *»*" 

PSIR0  =  PSID0*DEGTOR 
PHIR0  =  PHID0*DEGTOR 

PROCEDURAL  (P0 , Q0, R0 , BETAR0 , ALPHR0 , THETR0, ETAR0 
, GAMMAR=VT0 , PHIR0 , C) 

CALL  TRAP  (VT0 , PH1R0, G, GAMMAR , P0 , Q0 , R0, BETAR0 
, ALPHR0 . THETR0, ETAR0) 

END  5”  OF  TRAP  PROCEDURAL  " 

”***  CALCULATE  INERTIAL  CONSTANTS  ***" 


C0 

( (IXX*IZZ) - ( IXZ • IXZ) ) 

Cl 

= 

IZZ/C0 

C  2 

= 

IXZ/C0 

C3 

= 

C2* (IXX- IYY‘IZZ) 

C4 

= 

(  (IYY-IZZ) »C1)  -  ( I XZ* C2) 

C5 

= 

1  0/1 YY 

C6 

= 

C5* IXZ 

C7 

= 

CS*  (IZZ-IXX) 

C8 

= 

IXX/C0 

C9 

- 

( (IXX- IYY) *C8)  *  ( I XZ*C2) 
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DELETA=PULMAX*PULSE (ESTART, EREPET, EPULSE) 
DELPL  =THRMAX  *  PULSE (TSTART, TREPET, TPULSE) 
ETAR=ETAR0*DELETA 

CALL  CONTROLS (ETAR, DELPL, PLS) 

E TAD=ETAR » RADTOD 

END  $  "OF  CONTROLS  PROCEDURAL" 


**********  CALCULATE  QUATERNIONS;  NORMALIZE.  AND  THEN  ***»*»**»" 
"*“*“***  DETERMINE  THE  DIRECTION  COSINES  *»**«****•• 

PROCEDURAL (L1.L2,L3,M1.M2,M3,N1,N2, N3=TAU0 , TAU1 , TAU2, TAU3) 

"»**  NORMALIZE  QUATERNIONS  ***" 

TAUN  =  SQRT { (TAU0*TAU0) ♦ (TAU1*TAU1) ♦ (TAU2‘TAU2) ♦ (TAU3*TAU3) ) 

TAU0N  =  TAU0/TAUN 

TAU1N  =  TAU1/TAUN 

TAU2N  =  TAU2/TAUN 

TAU3N  =  TAU3/TAUN 

"***  CALCULATE  THE  DIRECTIONAL  COSINES  (LI  ->  N3)  ***" 

LI  =  (((TAU0N*TAU0N)*(TAU1N*TAU1N))*2.0)-1.0 

L2  =  ((TAU1N*TAU2N)+(TAU0N*TAU3N))*2.0 

L3  =  (  (TAU1N*TAU3N)  -  (TAU0N*TAU2N) )  *2 .0 

Ml  =  ((TAU1N*TAU2N) - (TAU0N*TAU3N) ) *2.0 

M2  =  ( ( (TAU0N*TAU0N) ♦ (TAU2N*TAU2N) ) *2 .0) -1 .0 
M3  =  ( (TAU2N*TAU3N) ♦ (TAU0N*TAU1N) ) *2 .0 

N1  =  ( (TAU1N*TAU3N) ♦ (TAU0N*TAU2N) ) *2 .0 

N2  =  ( (TAU2N*TAU3N) - (TAU0N*TAU1N) ) *2 .0 

N3  =  ( ( (TAU0N*TAU0N) + (TAU3N*TAU3N) ) *2 .0) -1 .0 

END  5 "OF  QUATERNION  PROCEDURAL" 


"**»  CALCULATE  ATMOSPHERIC  CONDITIONS  FOR  CURRENT  ALTITUDE  ***" 
CALL  ATMOS  (-Z.RHO) 


«  *  *  * 
"*  *  « 


CALCULATE  LONGITUDINAL  AND  LATERAL  AERODYNAMIC  FORCES 
IN  STABILITY  AXES  AND  MOMENTS  IN  BODY  AXES 


*  *  *« 
Hi" 


PROCEDURAL  (XA , YA , ZA , LA , MA , NA=VT , ALPHAR , BE  TAR , P , Q , R , PHI R , THETAR , DALPHR ) 

CALL  AERO  (XA. YA, ZA, LA, MA, NA, VT , ALPHAR , BETAR, P , Q, R, 

PHIR , THETAR , DALPHR) 

END  $"OF  AERO  PROCEDURAL" 


PROCEDURAL  (XP , YP , ZP , LP , MP , NP=VT , ALPHAR , BETAR , P , Q , R , PHIR , THETAR) 

CALL  PROP  (XP, YP,ZP,LP,MP,NP,VT, ALPHAR, BETAR, P,Q,R, 

PHIR, THETAR) 


END  S"OF  PROP  PROCEDURAL 


■T  -  ’  k  Ml 


SA  =  SIN (ALPHAS) 
CA  =  COS (ALPHAS) 
SB  =  SIN (BETAS) 
CB  =  COS (BETAS) 


”***  CALCULATE  GRAVITY  FORCE  COMPONENTS  ->  STABILITY  AXES  *»*" 

FGX  =  L3*MASS*G*CA  ♦  N3*MASS*G*SA 
FGY  =  M3‘MASS*G 

FCZ  =-L3*MASS*G*SA  ♦  N3»MASS*G*CA 


"***  CALCULATE  THE  TOTAL  FORCES  (FXS , FYS. FZS)  •**” 
"***  ->  STABILITY  AXES  ***" 

FXS  =  XA  ♦  XP *CA  ♦  ZP*SA  ♦  FGX 
FYS  =  YA  ♦  YP  ♦  FGY 
FZS  =  ZA  -  XP*SA  ♦  ZP»CA  ♦  FGZ 


"***  CALCULATE  THE  TOTAL  MOMENTS  (L, M, N)  ->  BODY  AXES  **•” 

L  =  (LA*CA)  -  (NA*SA)  ♦  LP 
M  =  MA  +  MP 
N  =  (LA*SA)  ♦  (NA*CA)  ♦  NP 


"**»  CALCULATE  THE  TOTAL  FORCES  (FXW, FYW, FZW)  ***" 
"***  ->  AIR-PATH  AXES  ***" 

FXW  =  FXS*CB  ♦  FYS*SB 
FYW  =-FXS*SB  ♦  FYS*CB 
FZW  =  FZS 


"***  CALCULATE  FORCES  IN  BODY  AXES  ***" 

FXB  =  XA*CA  -  ZA*  SA  ♦  XP 

FYB  =  YA  ♦  YP 

FZB  =  XA* SA  +  ZA*CA  ♦  ZP 


"***  CALCULATE  LINEAR  ACCELERATIONS  »**" 

AX  =  FXB/(MASS*G) 

AY  =  FYB/ (MASS*G) 

AZ  =  FZB/ (MASS*G) 


"***  CALCULATE  NORMAL  ACCELERATION  ***" 
AN  =  -AZ 


”***  RESOLVE  BODY  AXIS  ANCULAR  RATES  INTO  STABILITY  AXES  ***" 

PS  =  P*CA  -  R*SA 
RS  =  - P  *  SA  -  R  *CA 


******  CALCULATE  LINEAR  DERIVATIVES  ->  AIR  PATH  AXES 
***  DECLARE  DAL P HR  AS  AN  IMPLICIT  VARIABLE  ***" 


DAL P HR  =  IMPL  (Q,0  00 1.3(3. EE, 


(Q*CB-PS‘SB*FZW/MASS/VT)  /CB.  0 .  0001) 


PROCEDURAL  (=EF) 

CALL  IMP (EF) 

END  $"OF  DALPHR  PROCEDURAL" 


DBETAR  =  FYW/MASS/VT-RS 
DVT  =  FXW/MASS 


"***  CALCULATE  ANGULAR  ACCELERATIONS  ->  BODY  AXES  ***" 

PDOT  =  (L*C1) ♦ (N*C2) ♦ ( ( (P*C3) ♦ (R*C4) ) *Q) 

QDOT  =  (M*C5) ♦  {  ( (R*R) - (P*P) ) *C6) ♦ (R*P*C7) 

RDOT  =  (N*CB) + (L*C2) ♦  (  ( (P*C9) - (R*C3) ) *Q) 


DETERMINE  QUATERNION  RATES  AND  EULER  ANGLES. 


"***  CALCULATE  THE  QUATERNION  RATES  ***” 

TAU0DT=-  (  (  (TAU1N*P)  ♦  (TAU2N*Q)  +  (TAU3N*R) ) *0.5) 
TAU1DT=  ( (TAU0N*P)  -  (TAU3N*Q)  ♦  (TAU2N*R) )  *0 . 5 
TAU2DT=  ( (TAU3N*P) + (TAU0N*Q) - (TAU1N*R) )  *0.5 
TAU3DT=- ( ( (TAU2N*P) - (TAU1N*Q) - (TAU0N*R) ) *0.5) 


"***  CALCULATE  THE  NEW  EULER  .LWGLES  ***" 
PROCEDURAL (PSID. THETAD, PHID=L1 , L2 , L3 , M3 , N3) 


"»**  CALCULATE  HEADING  (OR  YAW)  ANGLE  (PSI)  ***" 
"***  IN  RANGE  0  (NORTH)  TO  360  DEG.  . 

PSIR  =  ATAN (L2/L1) 

IF (PSIR.LT.0.0)  PSIR  =  PSIR+PI * 2 
PSID  =  PSIR*RADTOD 


"***  CALCULATE  ANGLE  OF  PITCH  (THETA)  ***" 
"***  IN  RANGE  ♦  /-  90  DEGREES  ***" 

THETAR=  ATAN (-L3/SQRT ( (L1‘L1) ♦ (L2*L2) ) ) 
THE TAD =  THE TAR  *  RADTOD 


CALCULATE  BANK  (OR  ROLL)  ANGLE  (PHI)  IN  RANGE 
.  ♦/- 180  DEGREES, WHERE  0  DEG.  INDICATES  WINCS  LEVEL  ***" 

PHIR  =  ATAN (M3/N3) 

PHID  =  PHIR*RADTOD 

END  $"OF  EULER  PROCEDURAL" 


TRAJECTORY 


"•••  RESOLVE  FLICHT  PATH  AXES  VELOCITIES  INTO  BODY  *** 

"*•*  COMPONENTS  FOR  CALCULATION  OF  DISTANCES  IN  EARTH  AXES  *** 


UB  =  VT*CA*CB 
VB  =  VT*SB 
WB  =  VT*SA*CB 


"***  CALCULATE  EARTH  AXES  VELOCITIES  FROM  BODY  AXES  ***" 
"***  VELOCITIES  USING  DIRECTION  COSINES  AND  WIND  VELOCITIES  ***” 

VNE  =  L1*UB+M1*VB+N1‘WB-VWN 
VEE  =  L2*UB+M2*VB+N2*WB-VWE 
VDE  =  L3*UB+M3*VB*N3*WB-VWD 


PROCEDURAL (VGRKT , GAMMAD , CHID=VNE , VEE , VDE) 


"***  CALCULATE  THE  VELOCITY  OVER  THE  GROUND,  ***" 
"***  VGR  (IE.  GROUND  SPEED)  ***" 

VGR  =  SQRT((VNE*VNE) ♦ (VEE*VEE) ) 

VGRKT  =  VGR*MPSTOK 


"***  CALCULATE  FLIGHT-PATH  ANGLE,  GAMMA  IN  RANGE  +/-9(3  DEG  ***" 

GAMMAR  =  THETAR-ALPHAR 
GAMMAD  =  GAMMAR  *  RADTOD 


"***  CALCULATE  ANGLE  OF  CLIMB,  LAMBDA  IN  RANGE  */-9(3  DEG  ***" 

LAMDAR  =  ATAN ( - VDE /VGR ) 

LAMDAD  =  LAMDAR  *  RADTOD 


"***  CALCULATE  ANGLE  OF  TRACK,  CHI,  ***" 
"***  IN  RANGE  0 (NORTH)  TO  360  DEG.  ***" 

CHIR  =  ATAN (VEE/VNE) 

CHID  =  CHIR* RADTOD 

END  $"OF  TRAJECTORY  PROCEDURAL" 


"***  CALCULATE  THE  RANGE  (NOTE:  IN  METRES)  ***" 
RANGE  =  SQRT((X*X)*(Y‘Y)) 


INTEGRATION  OF  SYSTEM  STATE  EQUATIONS 


ALPHAR 

INTVC 

BE  TAR 

= 

INTVC 

VT 

= 

INTVC 

P 

= 

INTVC 

Q 

= 

INTVC 

R 

= 

INTVC 

(DALPHR . ALPHR0) 
(DBETAR, BETAR0) 
(DVT. VT0) 

(PDOT.P0) 

(QDOT.Q0) 

(RDOT.R0) 


$ "ROLL  RATE” 

$ "PITCH  RATE" 
$"YAW  RATE" 


TAU0  =  INTVC (TAU0DT, TAU00) 

TAU 1  =  INTVC (TAU IDT. TAU 10) 

TAU2  =  INTVC (TAU 2DT, TAU 20) 

TAU 3  =  INTVC (TAU 3DT, TAU 30) 

Z  =  INTEG (VDE . Z0) 

X  =  INTEG (VEE, X0) 


$"QUATERNION  TERMS" 


$ " Z - POS I T I ON AL  CO-ORDINATE" 
S "X  -  POSITIONAL  CO-ORDINATE" 


=  INTEC  (VNE ,  Y?) 


END  $"0F  DERIVATIVE" 


VE  =  VT*SQRT (RHO/RHO?) 

VEK  =  VE ‘MPSTOK 

ALPHAD  »  ALPHAR ‘RADTOD 

ALT  =-Z/FTTOM 


$"Y-POSITIONAL  CO-ORDINATE" 


$"EQUIVALENT  AIRSPEED  (M/S) " 
$ "EQUIVALENT  AIRSPEED  (KTS) " 
$ "ANGLE  OF  ATTACK  (DEG)" 
$ "ALTITUDE  (FT.)" 


EXPRESS  TERMINATION  CONDITION 


TERMT (TIME .GE . TSTOP) 


END  $"OF  DYNAMIC" 


"***  NOTE:  THIS  LISTING  SHOULD  BE  USED  IN  CONJUNCTION  WITH  ***" 

"***  THE  ADVANCED  CONTINUOUS  SIMULATION  LANGUAGE  (ACSL)  USER  ***" 
GUIDE/REFERENCE  MANUAL.  ***" 


END  $"OF  PROGRAM" 


rt 


SUBROUTINE  QUATNS (PS IR , THETAR , PHIR . TAU? , TAU1 , TAU2 , TAU3) 

C  "***  CALCULATES  NORMALIZED  QUATERNIONS  *»*" 

SPSI  =  SIN(PSIR*?5) 

CPSI  =  COS(PSIR*?.5) 

STHETA  =  SIN (THE TAR *0 . 5) 

C  THETA  =  COS(THETAR»?.5) 

SPHI  =  SIN  (PHIR*?. 5) 

CPHI  =  COS (PHIR*?. 5) 


C  "***  CALCULATE  THE  QUATERNIONS  FROM  THE  EULER  ANCLE  TERMS  ***" 

TAU?  =  (CPHI »CTHETA*CPSI) ♦ (SPHI *STHETA*SPSI) 

TAU1  =  (SPHI *CTHETA*CPSI) - (CPHI*STHETA*SPSI) 

TAU 2  =  (CPHI *STHETA*CPSI) ♦ (SPHI »CTHETA*SPSI) 

TAU 3  =  (CPHI »CTHETA*SPSI) - (SPHI *STHETA*CPSI) 


C  "***  NORMALIZE  INITIAL  QUATERNIONS  *»*" 

TAUN  =  SQRT( (TAU?*TAU?) ♦ (TAU1*TAU1) ♦ (TAU2*TAU2) 
♦ (TAU3*TAU3) ) 

TAU?  =  TAU?/TAUN 
TAU1  =  TAU 1/ TAUN 
TAU 2  =  TAU 2 /TAUN 
TAU 3  =  TAU3/TAUN 


RETURN 


SUBROUTINE  EVAL  (XXXX, XXXRES, MAXR) 


"**«  THIS  SUBROUTINE  IS  CALLED  FROM  THE  TRIMMING  SUBROUTINE 

TO  GAIN  ACCESS  TO  THE  DERIVATIVE  SECTION.  **•" 


DIMENSION  XXXX (3) , XXXRES (3) 

THE  H$"  CONTROL  CHARACTER  MAKES  MAIN  PROGRAM  ***" 
"***  VARIABLES  AVAILABLE  TO  THIS  SUBROUTINE.  ***" 


INTECER  I 
1=1 


.VV  >>«>. 
AV>'>‘V 


.  V  V  -  V'* 

5SSS& 


V. 


V  V 

WW1 


.  J  V  VJ 


THE TAR  =  (XXXX(l) ) 

ALPHAR  =  (XXXX (2)) 

ETAR0  =  (XXXX (3)) 

CALL  QUATNS  (PSIR0.THETAR , PHIRJ8,  TAU0. TAU1 .  TAU2,  TAU3) 

CALL  ZZDERV(I) 

XXXRES (1)  =  (DVT) 

XXXRES (2)  =  (DALPHR) 

XXXRES (3)  =  (QDOT) 

RETURN 

END 
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SUBROUTINE  IMP(EF) 


"*»»**  NOTIFIES  FAILURE  OF  IMPLICIT  ITERATION  ROUTINE  *«***" 


IF  (EF  .EQ.  0.JJ)  GOTO  89 
WRITE  (6.99) 

FORMAT (//27HIMPLICIT  PROCEDURAL  FAILED  //) 


-  -'A 


RETURN 
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FTSUBS.F 


SUBROUTINE  CAFCDI (PHIN , HNM , VN , IXX , I YY , I ZZ . IXZ . MASS) 

IMPLICIT  REAL (A-Z) 

INCLUDE  'FTPAR.F' 

INTECER  I.J 


C  ROUTINE  FOR  THE  INPUT  OF  CONFIGURATION  AND 

C  FLIGHT  CONDITION  DATA. 

OPEN  (UNIT=1, FILE= ' FTSUD2 .IN’, STATUS= ' OLD ' ) 

OPEN  (UNIT=7 , FILE= ' FTGRAF . VPP ' , STATUS= ' OLD ' ) 


C  ***  READ  IN  CONFIGURATION  AND  FLIGHT  CONDITION  DATA  *** 


READ  ( 1 .*) TTOT , TSAMP , NH . HN , DELH , NV , VN , DELV . NPHI , PHI N 

READ ( 1 , * ) DELPHI , NPL ,  PLp , DELPL , NRPM , RPM , DELRPM , WEGHT , NXCG, XCGP 

READ (1 ,  *) DELXCG, ZCG, (PEFF(I) ,1=1,7) , XT 

READ (1 , *) ZT, XTH, ZTH, THSET, MAXEP , NBLAD, PDIA, WSET,  TPSET,  ETAG 

READ (1 , * ) ACMASS , ACIXX . ACI YY . ACI ZZ , ACIXZ , SW. CW. BW, ST , CT 

READ ( 1 , * ) CETA, XP , ZP , CWP , BFW, BTAIL, CL0 , CLAL , CDp , CDAL 

READ ( 1 . ‘ ) CMP . CMAL , EPSP , EPSAL , QTOQ , Ap , A1 , A2 , A3 . Bp 

READ (1, *)  Bl, B2 , B3, CDPT , CDLT . CMTP , CMQW , CLP . CLXI , KWI NG 

READ ( 1 . * ) KFUSE , NPSF I , NPSFR , MPSF I , MPSFR , CDB , GTR , BLP , WTR , TAPERF 

READ ( 1 , * ) XWC , ZWC . SpS , XQARTC 

READ (7,*) (FPROP (1,1) , FPROP (I , 2) . 1=1 , 2p) 

READ (7 ,  *) ((TOP (I.J) , J=l, 33) ,I=l,2p) 

READ (7,*) ((ETP(I.J) ,J=l,6p) ,1=1,21) 

111  READ (7, *) (CYPROP (1,1) ,CYPROP(I,2) ,CYPROP(I,3) ,CYPR0P(I,4) ,I=l,2p) 
READ (7 , *) ( (DELEPS(I, J) ,J=l,ip) ,1=1,8) 


MASS=ACMASS 

IXX=ACIXX 

IYY=ACIYY 

IZZ=ACIZZ 

IXZ=ACIXZ 

HNM=HN 

CLOSE (UNIT=1) 
CLOSE (UNIT=7) 


RETURN 

END 


SUBROUTINE  PARINC  (VECIN , CCPOS , PLS) 

C  "  SUBROUTINE  TO  ALLOW  RUNTIME  PARAMETER  MODIFICATIONS  " 

IMPLICIT  REAL (A-Z) 

DIMENSION  VECIN (6) 

INCLUDE  'FTPAR.F' 
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ALTp=  VECIN(l) 
RPM  =  RPM-VECIN(2) 
XCGP=XCGP ♦ VECI N (3) 
PLp  =PLp  +VECIN (4) 
PAR5=PAR5*VECIN(5) 
PAR6=PAR6«-VECIN (6) 


XCG=XCGP*CW/ipp . ♦XQARTC-CW/4 . 

CGPOS=XCCP 

PLS=PLP 

RETURN 

END 


SUBROUTINE  TRAP  (VTp , PHIRp, G,  GAMMAR ,  P p ,  Qp ,  Rp ,  BETARP 
,  ALPHRp ,  THETR0 ,  ETARp) 

"  GIVES  INITIAL  APPROXIMATION  FOR  TRIM  CONDITION 

IMPLICIT  REAL (A-Z) 

INCLUDE  ' FTPAR . F ' 


PI  =  3.141593 
ETAPP=  2.2 

KK  =41 .p* (XCGP/ipp . -p.4) 

QD  =  ( . 5*1 . 2256* VTp* VTP) 

PP  =p. 

QP  =G/VTP*TAN (PHIRp) *SIN (PHIRP) 
RP  =G/VTP*SIN (PHIRP) 

AN  =  1/COS (PHIRP) 

BETARp  =  p. 


ALPHAW  =  ACMASS*G/(CLAL*QD*SW) 
ALPHRp  =  ALPHAW  -  WSET 
CLWB  =  CLP  ♦  (CLAL* ALPHAW)  +  CLF 
CDWB  =  CDP  ♦  (CDAL*CLWB**2)  ♦  CDF 

DRAG  =  CDWB  *  QD  *  SW 

GAMMAR  =  ASIN (-DRAG/ (ACMASS*G) ) 
THETRp  =  GAMMAR  *  ALPHRP 

ETAP  =  (ETApp»KK*CLWB) »PI/18p. 
ETARP  =  ETAp 


RETURN 

END 


SUBROUTINE  ATMOS  (H.RO) 
IMPLICIT  REAL (A-Z) 


V  w  'w-  *V  V.  -A' l".  1PJ  7;  jr. .-; .".  7-.  -• 


C0MM0N/ARGS2/HITE, DENS 
HITE=H 

DENS=1. 2256* (1-0. 02256*HITE/1000.0)  “4.2561 
RO  =DENS 


RETURN 

END 


SUBROUTINE  AERO  (XA, YA, ZA, LA, MA, NA.VT. ALPHAR. BETAR 
1  ,PP.QQ,RR,PHIR.THETAR.DALPHR) 

C 

C  MODIFIED  BY  RODD  PERRIN  OCTOBER  1985 
C 


IMPLICIT  REAL (A-Z) 

INCLUDE  'FTPAR.F' 

DIMENSION  X (20) 

XDAL  =  DALPHR 

X  (V)  =  VT 
X (ALPHA)  =  ALPHAR 
X(Q)  =  QQ 
X(P)  =  PP 
X(H)  =  HITE 
X  (PHI)  =  PHIR 
X  (THETA)  =  THE  TAR 

CALL  ATMOS (HI TE.RHO) 

QD=.5*RHO*VT**2 

CALL  THRUST (X. TOP, ETP) 

CALL  PROPEL (X) 

CALL  FLAPS (X) 

CALL  WB  (X) 

CALL  TAIL(X) 

CDW^CDWB - CD  B 

CALL  RESOLV(X (ALPHA) , CLWB, CDW, XP, ZP . CMWBF) 
ALTF =ALPHAT - TP  B 

CALL  RESOLV (ALTF , CLT , CDT , XT , ZT , CMTF ) 


CTHXW=CTPROP  *COS (ALTH) -CLPROP*SIN (ALTH) 
CTHZW= -CLPROP  *COS (ALTH) -CTPROP»SIN (ALTH) 

CTHX=CTPROP  *  COS (THSET) -CLPROP‘SIN (THSET) 
CTHZ=-CLPROP*COS (THSET) -CTPROP*SIN (THSET) 

ALTFS=X (ALPHA) -ALTF 

CTXW= -CLT*SIN (ALTFS) -CDT*COS (ALTFS) 

CTZW=-CLT*COS (ALTFS) *CDT‘SIN (ALTFS) 

CMTHF=CTHX‘ (ZTH-ZCG) ♦CTHZ* (XTH-XCG) 
CLWBT=CLWB- ST* (CTZW) /SW 

CALL  RESOLV (X (ALPHA) ,0.0, CDB , XP , ZTH , CMBF ) 
CMWBTH=CMWB* (CMWBF  +CMTHF  «CMBF)  /CW 
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CALCULATE  AERODYNAMIC  FORCES 
XA  =  QD#  (GW#  (CTIIXW  CDWP)  -  CT^CTXW) 


ZA  =  QD* (SW» (CTHZW-CLWBT) ) 


=  9-9 

-  QD* (SW*CW*CMWBTH* (ST* (CT*CMT*CMTF) ) ♦ (SW*CW*CMQW*Q) ) 

=  99 


RETURN 

END 


SUBROUTINE  PROP  (XP ,  YP ,  ZP ,  LP ,  MP ,  NP ,  VT , ALPHAR , BE TAR 
. PP , QQ , RR , PHIR . THETAR) 


RETURN 

END 


SUBROUTINE  CONTROLS (ETAR , DELPL , PLS) 

INCLUDE  'FTPAR.F' 

ETA  =  ETAR 
PL  =  PLJJ+DELPL 
PLS  =  PL 


RETURN 


FTSUBS2.F 


SUBROUTINE  THRUST(X) 

C  ! CALCULATES  COMPONENTS  OF  PROPELLOR  THRUST 

IMPLICIT  REAL  (A-Z) 

INCLUDE  'FTPAR.F' 

INTEGER  I.J 
DIMENSION  X(20) 

1000  FORMAT ('  POWER  TOO  LOW') 

1001  FORMAT ('  POWER  TOO  HIGH') 

4  EP=PL*MAXEP 

ETAP=0 
TCC=0 .0 
TTHST =0.0 
IF (RPM==0) CO  TO  1 
IF  (PL==0) GO  TO  1 

IF (X (V) >55 .0) GO  TO  5 

POD2=S0S  *  EP/745 . 7/ ( (PDIA/0 . 3048)  *  *  2) 

PIND=3 . 142  *RPM*PDIA/0 . 3048/60 . 0 

DO  3  J=l, 6 
DO  2  1=1,5 

2  CALL  INTERP  (TOP (1,1) ,TOP(l,  (J-l) *5* I +3) , 20,POD2,TOP (1,34)  , ERROR) 

3  CALL  INTERP  (TOP  (1,2) , TOP (1, 34) , 5 , PIND, TOP (J, 35) , ERROR) 

CALL  INTERP (TOP (1,3) , TOP (1 , 35) , 6,X (V) ,TDP, ERROR) 

IF (TDP<0) GO  TO  9 
IF (TDP>10) GO  TO  10 
TTHST=TDP*EP* 4 .4497/745 . 7 
TCC=TTHST/(2.0*QD*PDIA**2) 

ETAP=TTHST*X (V) /EP 
RETURN 

5  S0SCP=S0S*EP/DENS/  ( (RPM/60) *  *  3) / ( (PDIA) *  *5) 

IF (S0SCP<0.025)GO  TO  9 

MACH=X (V) /340 . 29 

PINDOA=3 . 142 ‘RPM‘PDIA/60 . 0/340 . 29 

DO  8  J=l, 5 
DO  7  1=1,11 

7  CALL  INTERP (ETP (1,1) ,ETP(1,  (J-l) *11  +  1-3) , 21, S0SCP, ETP  (1 , 59) 

1  , ERROR) 

8  CALL  INTERP  (ETP  (1, 2)  , ETP (1, 59)  , 11 , PINDOA, ETP ( J, 60)  , ERROR) 

CALL  INTERP  (ETP  (1, 3) , ETP (1,60) , 5, MACH, ETAP, ERROR) 

IF (ETAP>1  0)GO  TO  9 
IF (ETAP<0 . 0) GO  TO  11 
1  TTHST=ETAP *  EP/X (V) 

TCC=TTHST/(2 .0*QD*PDIA* ‘ 2) 

RETURN 

9  IF  (RPM<600) GO  TO  12 
RPM=RPM- 200 

GO  TO  4 

12  TTHST= -  1  . 0 

TYPE  1000 
RETURN 

10  IF (RPM>2600)GO  TO  13 


13 


14 


RPM=RPM*200 
GO  TO  4 
ETAP=-1 .0 
TYPE  1001 
RETURN 

IF (RPM<600) GO  TO  14 

RPM=RPM-200 

GO  TO  4 

TCC= -1.0 

TYPE  1001 

RETURN 

END 
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SUBROUTINE  WB (X) 

C  ! CALCULATES  WING- BODY  AERODYNAMICS 

IMPLICIT  REAL  (A-Z) 

INCLUDE  ' FTPAR . F ' 

DIMENSION  X (20) 

CLWB=CL0* (CLAL*ALPHAW) *CLF 
CDWB=CD0+ (CDAL* (CLWB*  *  2)  )  +CDF 
CMWB=CM0+ (CMAL * ALPHAW)  *CMF 

CLSS=0 

CMSS=0 

QS0QIN=TCC*8 .0/3 . 142 
TEMP=DUWDAL 

CALL  SLIPST (X, 2) 

IF (ZSS> (PDIA/2 . 0) ) GO  TO  2 

BWP=2*SQRT ( (PDIA/2 .0) *  *2~ZSS*  *2) -BFW/2 

IF  (PEFF  (2) ==0)GO  TO  1 

CMSS= (CM0+CMF ) *CWP  *  »  2  *  BWP  «QSOQIN/CW/SW 
1  IF  (PEFF (3) ==0) GO  TO  2 

DELAW=- (DPEDAL*ALTH) / (1+TEMP) 

AIL=BWP/CWP 

SIL=BWP*CWP 

BL=0 . 0B 

ARL=BW/CW 

IF  (ARL>10.0) BL=0.0 

K2=3 . 0* (AIL-1.0) +BL* (10.0-ARL) *AIL 

K3=1.96*SQRT(TCC) 

K4=  (K2  +  3 .0) *0.1*K3 

K1=0 . 2*K4+1 . 181-0 . 489*SQRT (TCC*0 . 1) 

CLSS=K1*SIL* ( (1+QSOQIN) *CLAL*DELAW*QSOQIN*CLWB) /SW 


*  * ' 


2  CLWB=CLWB+CLSS 
CMWB=CMWB ♦ CMS  S 

CDWB=CDWB* (1 . 0+CLSS/CLWB) 

XZSS (18) =ZWC 
XZSS (19) =ZT 
XZSS (20) =PDIA 

RETURN 

END 


SUBROUTINE  TAIL(X) 

C  (CALCULATES  TAIL  CONTRIBUTION 

IMPLICIT  REAL (A-Z) 

INCLUDE  ' FTPAR. F' 

DIMENSION  X (20) 
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CALL  BEND (X) 


CALL  SLIPST (X, 3) 

IF  (ZHEFF>  (PDIA/2  .0)  )G0  TO  1 
SHI=2*CT»SQRT ( (PDIA/2 .0) * ‘2-ZHEFF*  *2) 

K 1=4. 08* SHI /ST 
K2=K1* (TCC*-0. 575) -0.3 

K3=0. 202-0. 03*ZHEFF/  (PDIA/2. 0) -0.0754* ( (ZHEFF/ (PDIA/2 .0) ) * *2) 
DELQT=K2*K3-0.1 
IF (TCC<0 . 1) DELQT=0 . 0 
GO  TO  2 

DELQT=0 
RT=QTOQ*DELQT 
RT1=RT- 1 

IF  (PEFF (5) ==0) RT*1 .0 

ALPHAQ= (XT-XCG) *X(Q)/(X(V) *SQRT (RT) ) 

ALPHAT=X (ALPHA) +TPB-EPS+ALPHAQ 
CLT=RT* (A0*A1*ALPHAT) +A2*ETA*A3*BETA 
CH=RT* (B0*B1*ALPHAT) ♦B2*ETA*B3*BETA 

IF  (PEFF (6) ==0) RT1=0 . 0 

CLT=CLT+RT1* (A2*ETA+A3*BETA) 

CH=CH+RT1* (B2*ETA+B3*BETA) 

CDT=CD0T*CDLT*CLT* *  2 

CLTTH=CLT-RT* (A2  *  E  TA  *A3  *  BE  TA) 

CDTTH=CD0T +CDLT  *  CLTTH* *  2 

CMT=CMT0 

PETA=QD*CETA*CH*ETAG 

XZSS  (22) =CT 

RETURN 

END 


SUBROUTINE  RESOLV (AL.CLl, CD1 , XI , Z1 , CM1) 
! TRANSLATES  FORCES  TO  THE  CG 
IMPLICIT  REAL (A-Z) 

INCLUDE  'FTPAR.F' 

CM1  =  -  (CD1 *COS  (AL)  -CH*SIN  (AL) )  *  (Zl-ZCG) 
1  - (CLl'COS (AL) +CD1*SIN (AL) ) « (Xl-XCG) 

RETURN 

END 


SUBROUTINE  DWASH (X) 

! CALCULATES  DOWNWASH  AND  TOTAL  HEAD  AT  THE  TAIL 
! AND  ALSO  INDUCED  INCIDENCE  AT  THE  TAIL  DUE  TO  PITCHING 
IMPLICIT  REAL (A-Z) 

INTEGER  I 
INCLUDE  ’FTPAR.F' 

DIMENSION  X ( 20) 

KA=CW/BW- 1 .0/(1 .04 (BW/CW) **1.7) 

KLAM- ( 10 . 0- 3 . 0*WTR) /7 . 0 
LH  =  XT - XWC 

HH= (XT  - XWC)  *TAN (WSET) *ZWCZT 
KH= (1 -HH/BW) / ( (2 ,0*LH/BW) **0.33) 


V,Vv  . 

ft  V*  »  * 


>v.*v 


O'.- 


.•  V  >  V  J 


■v'.V.'v' 


'  V.sV 


EPSAL=4 . 44* ( (KA*KLAM*KH) **1.19) 

EPS=57 . 3* (EPS0+EPSAL* (ALPHAW-XDAL* (XT-XCG) /X (V) ) ) 

IF(PEFF(4)==0)CO  TO  1 
DO  2  1=1,8 

CALL  INTERP (DELEPS (1, 1) , DELEPS (1 , 1*2) , 8, EPS.DELEPS  (1 , 11) .ERROR) 
CALL  INTERP  (DELEPS  (1,2) , DELEPS  (1. 11) , 8, TCC, DEP, ERROR) 

ZHT=ZTH» (XT-XTH) *TAN (THSET) -ZT 
KFAC=0 . 85-0 . 25*ABS (ZHT) /PDIA 
DEPS=KFAC*DEP 
EPS=EPS+DEPS 

EPS=EPSF+ (EPS)  /57 . 3 
QTOQ=QTOQ 

ALPHAQ= (XT-XCG) *X (Q) /X (V) 

XZSS  (21) =CW 

RETURN 

END 


SUBROUTINE  BEND (X) 

! CALCULATES  TAILPLANE  ANGLE  DUE  TO  FUSELAGE  BENDING 
IMPLICIT  REAL (A-Z) 

INCLUDE  'FTPAR.F ’ 

DIMENSION  X (20) 

LT=QD*QTOQ*ST*CLT 

TPB=KFUSE*LT+TPSET 

RETURN 

END 


SUBROUTINE  TORSON (X) 

! CALCULATES  EFFECTIVE  WING  INCIDENCE 
IMPLICIT  REAL (A-Z) 

INCLUDE  'FTPAR.F' 

DIMENSION  X (20) 

LW=QD*SW*CLWB 

ALPHAW=X (ALPHA) +WSET-KWING*LW 

RETURN 

END 


SUBROUTINE  PROPEL (X) 

•CALCULATES  PROPELLER  FORCES 

IMPLICIT  REAL (A-Z) 

INCLUDE  'FTPAR.F' 

DIMENSION  X (20) 

INTEGER  IB 
CALL  TORSON  (X) 

CALL  SLIPST  (X , 1) 

ALTH=THSET«X  (ALPHA)  *DUWDAL ‘ALPHAW 
IF (RPM==0) GO  TO  1 


DPEDAL=0 


nnnnnnonn 


CLPROP=0 

IF (PEFF (1)==0)GO  TO  3 
IB=IFIX (NBLAD) 

CALL  INTERP {CYPROP (1, I) , CYPROP ( 1 , I B) , 20, BL0, CYPSI0. ERROR) 

CALL  INTERP (FPROP (1,1) ,FPROP(l,2) , 20, TCC.FPR, ERROR) 
CLPROP=3 . 142  *PDIA*  * 2  *CYPSI0«FPR»ALTH/SW/4 . 0 
C1=SQRT ( (SQRT (50 . 6*TCC+28 . 623) -5.35) /25 . 3) 

C2=0 . 27*EXP (-0 . 127*TCC) 

DPEDAL=C1+C2*CYPSI0 
3  CTPROP=  2  *  TCC  *PD I A*  *  2 /SW 

GO  TO  2 

1  DPEDAL=0 
CLPROP=0 
CTPROP=0 

2  CONTINUE 

XZSS (17) =ZTH 

RETURN 

END 


SUBROUTINE  INTERP (XLIST, YLIST, N, X, Y, ERROR) 

USED  FOR  INTERPOLATING  INTO  A  ONE  DIMENSIONAL  ARRAY 
X:  INDEPENDENT  VARIABLE 

XLIST:  LIST  OF  INDEPENDENT  VARIABLE  BREAK  POINTS  IN 

ASCENDING  ORDER 

YLIST:  LIST  OF  DEPENDENT  VARIABLE  BREAK  POINTS 

N:  NUMBER  OF  BREAK  POINTS 

Y:  INTERPOLATED  VALUE  OF  DEPENDENT  VARIABLE 

ERROR:  .TRUE.  IF  EXTRAPOLARION  WAS  NEEDED 

DIMENSION  XLIST (100) , YLIST (100) 

IF (X.GE .XLIST (1) ) GO  TO  5 
ERROR= . TRUE . 

1=1 

GO  TO  20 
5  IF (X.LE .XLIST (N) )GO  TO  7 

ERROR= . TRUE . 

I=N-1 
GO  TO  20 

7  ERROR= . FALSE . 

DO  10  1=1, N-l 

IF(X.GE.XLIST(I) .AND. X.LE. XLI ST (I ♦ 1) ) GO  TO  20 
10  CONTINUE 

I=N- 1 
20  J=I 

XLB=XLIST(J) 

XUB=XLIST (J+l) 

RX=XUB-XLB 
Wl=  (XUB-X) /RX 
W2= (X-XLB) /RX 

Y=W1*YLIST(J) *W2  * YLIST ( J* 1) 

RETURN 

END 


SUBROUTINE  SLIPST(X.I) 

C  (CALCULATES  LOCUS  OF  PROPELLER  SLIPSTREAM 


IMPLICIT  REAL (A- Z) 
INCLUDE  ' FTPAR . F ' 


DIMENSION  X (20) 
INTEGER  I 


GO  TO (10, 20, 30) , I 
XZSS (3) =XTH 

DUWDAL=0 . 28/ (XWC-XTH) /CW 
XZSS  (1) =2*XZSS (3) 

ALUW=X (ALPHA) ♦ DUWDAL  *  ALPHAW 
IF  (PEFF  (7) ) 12 , 12 , 11 
ALUW=X (ALPHA) + DUWDAL* (CLWB/CLAL) 
XZSS  (2) =ABS (XZSS (3) ) *TAN (ALUW) ♦  ZTH 
XZSS (4) =ZTH 

RETURN 


XZSS  (5) =XZSS (3) /2 
ALUW=ALUW- DPEDAL * ALTH 

XZSS (6) =  -ABS (XZSS (5) ) *TAN (ALUW) +XZSS (4) 

DUWDAL =0 . 28/ (XWC- (XTH/2) ) /CW 
IF  (PEFF  (7))21,21,22 
ALUW=ALUW+DUWDAL* ALPHAW 
GO  TO  23 

ALUW-ALUW+ DUWDAL* (CLWB/CLAL) 

XZSS (7) =XWC 

XZSS (8) =-ABS (XZSS (5) -XZSS (7) ) *TAN (ALUW) ♦XZSS (6) 
ZSS=ABS (ZWC-XZSS (8)  ) 

ALUW=-WSET 

XZSS (9) =OCWC+0 . 75*CW 

XZSS (10) =-0 . 75*CW*TAN (ALUW) ♦XZSS (8) +ZFLAP 
RETURN 


KA=CW/BW- 1/ (1+ (BW/CW) **1.7) 

KLAM= (10.0-3 .0*WTR) /7 . 0 
WTE=XWC*0 . 75*CW 

XZSS (11) =WTE* (XT-WTE)  /2 
LH=XZSS (11) -XWC 

HH= (XZSS (11) -XWC) *TAN (WSET) ♦ZWC-XZSS (10) 

1  - (XZSS (11) -WTE) *TAN (ALPHAW) 

KH=(1-HH/BW)/((2*LH/BW) *»0.33) 

DDWDAL=4 . 44* ( (KA*KLAM*KH) **1.19) 

XZSS (12) =  -ABS (XZSS (11) -WTE) *TAN (X (ALPHA) - (DDWDAL ‘ALPHAW) -EPSF) 
1  *XZSS(10) 

XZSS (13) =XT 
LH=XT-XWC 

HH=  (XT-XWC) *TAN (WSET) ♦ ZWC -XZSS (1 2) 

1  -  (XT- XZSS (11) ) ‘TAN (ALUW) 

KH=  (1-HH/BW)  /  (  (2*LH/BW)  “0 . 33) 

DDWDAL  =  4 . 44*  (  (KA‘KLAM‘KH)  “1.19) 

XZSS (14) =-ABS (XT-XZSS (11) ) ‘TAN (X (ALPHA) - (DDWDAL ‘ALP HAW) -EPSF) 

1  *  XZSS  (12) 

ZHEFF=ABS (ZT-XZSS (14)  ) 

XZSS (15) =XT* (XT-WTE) /2 
LH=XZSS (15) -XWC 

HH= (XZSS (15) -XWC) ‘TAN (WSET) ‘ZWC 
1  -XZSS (14) - (XZSS (15) -XT) ‘TAN (ALUW) 

KH  -  ( 1  -  HH/BW)  /  (  (2  •  LH/BW)  “0.33) 


DDWDAL=4 . 44* ( (KA*KLAM*KH) **1.19) 

XZSS (16) =-ABS (XZSS (15) -XT) *TAN (X (ALPHA) - ( DDWDAL  *  ALPHAW) -EPSF) 
1  ♦XZSS (14) 

RETURN 

END 

SUBROUTINE  FLAPS (X) 


! CALCULATES  CL, CM, CD. EPS, ZFLAP  FOR  A  SINGLE  SLOTTED 
!FLAP  FITTED  TO  THE  A1 .DEFLECTIONS  20  DEG  AND  40  DEG 


IMPLICIT  REAL (A-Z) 

INCLUDE  'FTPAR.F ' 

CLF=0.0 
CMF  =0 . 0 
CDF=0.0 
ZFLAP=0 .0 

IF (PEFF (7) -1) 15, 5, 10 

CLF=0 . 74 
CMF=-0 . 167 
CDF=0.0213 
ZFLAP=0 . 2*CW*0 . 364 
GO  TO  IS 

CLF=1 . 13 
CMF=-0 .2829 
CDF=0 . 0628 
ZFLAP=0 . 2*CW*0 . 839 

EPSF=  (CLF* (10 . 3+64 . 1* ( (ABS (ZWC-ZT) / (BW/2 .0) ) -0 . 53) *  *2) / 
1  (2.0»BW*5.0/SW))/57.3 
XFLAP=0. 2*CW 

RETURN 

END 
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PARAMETER  (V=l , ALPHA=2 , Q=3 . P=4 , H=5 , PHI=6 , THETA=7) 

INTEGER  V, ALPHA, Q,P ,H, PHI , THETA 

COMMON/ARGS/ACMASS , ACI XX , AC I YY , AC I Z  Z , AC IXZ , SW , CW , BW . WSE  T 
1  , ST, CT , CETA, TPSET , ETAG , XP , ZP , XT, ZT, XQARTC 

1  ,XTH,ZTH,ZSS,XCG,XCGP,ZCG,A0,A1,A2.A3,BJI,B1,B2,B3 

2  , CL0 , CD0 , CM0 , CLAL , CDAL , CMAL , MAXEP , THSET, ALTH 

3  ,  DUWDAL ,  ALUW,  EPS0,  EPSAL,  QTOQ,  KFUSE ,  KWING,  CD0T,  CDLT,  CMTJJ 

4  , CMQW, CLP.CLXI , TSAMP , V0, PHI0 

5  ,  PW, QW , RW , CTHX , CTHZ , CLWB , CDWB , CMWB , CLWBT 

6  , ALPHAW, TPB , ALPHAQ , EPS , DP E DAL , ALPHAT , CLT, CDT , CMT , CH 

7  , ETA0, ETA.PETA, BETA, XI 

8  , TTHST, CTHXW, CTHZW, CMWBF , CMTF , CLTTH, CDTTH, CMWBTH, CMTTH 

9  , NPSE I , NPSFR ,MPSFI,MPSFR,CDB, GTR , QD , RHO , XDV , XDAL 

1  , NBLAD , PDI A , RPM , CWP , BFW, CLPROP , CTPROP , ETAP , TCC , POP 

2  , WTR , TAPERF , BTAIL, DTAIL , BL0 , S0S , DEPS , CLF , CMF , CDF , EPSF 

3  , XFLAP , ZFLAP, DELEPS (8 , 11) , FPROP (20, 2) , CYPROP (20, 4) ,PEFF(7) 

4  , XZSS  (24)  ,XWC,ZWC.ZHEFF,RT,LT.PL,PL0,TOP(20,35)  ,ETP(21,60) 

COMMON/ARGS 2/HI TE , DENS 


LOGICAL  ERROR 


APPENDIX  3.  SET  UP  PROGRAM  LISTING  -  FTCHOO 


ETCHOO.F 


C  !  PROGRAM  FOR  PREPARING  DATA  FOR  INPUT  INTO  CAM’S  FLIGHT 

C  !  SIMULATION  PROGRAM. THE  DATA  IS  READ  IN  FROM  FTSUD  FILE 

C  !  AND  THEN  DISPLAYED  ON  A  VDU  ALONG  WITH  VARIABLE  NAMES 

C  !  CHANGES  CAN  BE  MADE  BY  APPLYING  THE  COMMANDS  INDICATED 


DIMENSION  A  (40)  ,B(54) 

DOUBLE  PRECISION  A1 (40) , B1 (54) 


DATA  AI/’TTOT',  'TSAMP',  ’NH\  'HN\  ' DELH'  ,  'NV\  'VN',  ' DELV '  , 
1  'NPHI ' . 'PHIN' , 'DELPHI ' , ' NPL ' , ' PLN ' , ' DELPL ' , 'NRPM' , 

1  'RPMN* , ' DELRPM ' . 'WEICHT' . 'NXCG' , 'XCG  %' , 'DELXCG' , 

1  ' ZCG' . 'PEFF(l) * , 'PEFF (2)  ' , ’PEFF(3) ' , ’PEFF(4)  ' , 'PEFF  (5)  ' , 
1  'PEFF  (6)  ' , 'PEFF (7)  * , 'XT', 'ZT\  'XTH' , ' ZTH* ,  'THSET' , 

1  'MAXEP ' , ' NBLAD ' , 'PDIA' , 'WSET' , 'TPSET' , 'ETAG'/ 


DATA  BI/ 'MASS ' , 'IXX', ' IYY' , 'IZZ', ' IXZ ' , 'SW'. 'CW' 

1  . 'BW' , 'ST' , 'CT' , 'CETA' , 'XP' , 'ZP' , 'CWP' 

1  , ' BFW' , ' BTAIL ' , 'CL0', 'CLAL'. 'CD0* , ' CDAL ' , 'CM0' 

1  . 'CMAL' , 'EPS0' , 'EPSAL' , 'QTOQ' , 'A0’ , 'Al' , 'A2' , 'A3' , 'B0' , 'Bl 
1  ,  'B2\  'B3\  '  CD0T ' ,  '  CDLT ' ,  '  CMT0 ' ,  'CMQW',  'CLP',  'CLXI\  'KWING 
1  , ' KFUSE ' , ' NPSFI ' , 'NPSFR' , 'MPSFI', 'MPSFR' , ' CDB ' , 'GTR', ' BL0 ' 
1  ' WTR ' , 'TAPERF', 'XWC', 'ZWC', ' S0S ' , ' XQARTC ' / 


OPEN  (UNI T=1 , F I LE= ' FTSUD2 . OUT ' , STATUS= ' OLD ' ) 

OPEN  (UNIT=2 , FILE= 'FTSUD20 .OUT ' , STATUS= 'OLD ' ) 

C  !  PRIMARY  DATA  READ  IN, STORED  IN  ARRAY  A 

READ (1 , *) (A (I) ,1=1,10) 

READ (1, *) (A (I ) ,1=11,20) 

READ ( 1 , * ) (A (I) ,1=21,30) 

READ(1, *) (A (I) ,1=31,40) 

C  !  SECONDARY  DATA  READ  IN, STORED  IN  ARRAY  B 

READ (1 , *) (B (I) , 1=1,10) 

READ (1, *)  (B  (I) , 1=11,20) 

READ (1 , *) (B (I) , 1=21,30) 

READ (1,  *) (B (I ) ,1=31,40) 

READ (1 , *)  (B  (I) , 1=41,50) 

READ (1, *)  (B  (I) ,1=51,54) 

C  !  FLAG  SETTING  -  0.0  FOR  ARRAY  A  .  1.0  FOR  B 

FLAG=0 . 0 

C  !  DISPLAY  ON  VDU  HEADINGS , THEN  NAMES  AND  VARIABLE  VALUES 

1502  PRINT  100 
PRINT  101 
PRINT  102 

IF (FLAC==1 .0)  GO  TO  1700 

DO  1000  1=1,20 
J=I ♦ 20 

PRINT  103, I ,A1 (I) ,A(I) . J, Al (J) ,A(J) 

1000  CONTINUE 


1000  CONTINUE 


1504  PRINT  104 

1500  ACCEPT  * . I ,  F 

IF (FLAG==1 .0)  CO  TO  1800 

A(I)=F 

GO  TO  1801 

1800  B ( I ) =F 

1801  CONTINUE 

IF (I  .NE.  0)  GO  TO  1500 

1501  PRINT  105 
ACCEPT  *.I 

GO  TO  (1502,1600,1503)  1*2 


1503  FLAG=1.0 

GO  TO  1502 

1700  DO  1001  1=1,27 
J=I +27 

PRINT  103 ,1,B1(I) , B  (I )  , J , B1  ( J)  , B ( J) 
1001  CONTINUE 

CO  TO  1504 


c  !  PRIMARY  DATA  WRITTEN  OUT  FROM  ARRAY  A 

1600  WRITE (2, *) (A(I) .1=1.10) 

WRITE (2, ‘) (A (I) ,1=11,20) 

WRITE  (2 , *)  (A (I) , 1=21,30) 

WRITE  (2,*)  (A (I )  ,1=31,40) 

c  !  SECONDARY  DATA  WRITTEN  OUT  FROM  ARRAY  B 

WRITE (2, *) (B(I) ,1=1,10) 

WRITE  (2 , *)  (B (I) ,1=11,20) 

WRITE (2 , *) (B (I) ,1=21,30) 

WRITE (2 , *) (B(I) .1=31,40) 

WRITE (2, *) (B (I) ,1=41,50) 

WRITE (2, ♦)  (B (I ) .1=51.54) 


100  FORMAT (/, /, 9X, '  '***  FTSIM  SETUP  DATA- TO  CHANGE  TYPE  I,F' 

1,  '  *  *  *  *  ' ) 

101  FORMAT (/, 17X, '  -  TERMINATE  CHANGES  WITH  0,0  -  ') 

102  FORMAT (/, 2 (’  #  (I) ' , IX, 'QUANTITY ' , 7X, 'VALUE (F) ' , 4X  )) 

103  FORMAT (2 (IX, 1 2 , 3X, A10 , F 14 . 7 , 2X) ) 

104  FORMAT (/,'  TO  CHANGE  VALUE  TYPE  I , NEW  F  ’) 

105  FORMAT (/, '  TYPE  -1,0,1  TO  REVIEW, TERMINATE  OR  VIEW  REMAINING' 
1 , '  DATA  ' ) 


CLOSE (UNIT=1) 
CLOSE (UNIT=2) 


FTCHOO 


““  FTSIM  SETUP  DATA-TO  CHANGE  TYPE  I.F  **** 


-  TERMINATE  CHANGES  WITH  0.0 


#(I) 

QUANTITY 

VALUE (F) 

#(I) 

QUANTITY 

VALUE (F) 

1 

TTOT 

21 

DELXCG 

.1135060 

2 

TSAMP 

. 1000000 

22 

ZCG 

.0591000 

3 

NH 

1 .0000000 

23 

PEFF(l) 

1 . 0000000 

4 

HN 

1.0000000 

24 

PEFF  (2) 

1.0000000 

5 

DELH 

.0000000 

25 

PEFF  (3) 

1 . 0000000 

6 

NV 

21.0000000 

26 

PEFF (4) 

1 . 0000000 

7 

VN 

60 . 0000000 

27 

PEFF (5) 

1 . 0000000 

8 

DELV 

3.0000000 

28 

PEFF  (6) 

1 . 0000000 

9 

NPHI 

1.0000000 

29 

PEFF (7) 

. 0000000 

10 

PHIN 

. 0000000 

30 

XT 

3.0000000 

11 

DELPHI 

100 . 0000000 

31 

ZT 

- .5500000 

12 

NPL 

. 0000000 

32 

XTH 

-3.6600001 

13 

PLN 

. 0000000 

33 

ZTH 

. 0000000 

14 

DELPL 

. 0000000 

34 

THSET 

. 0000000 

15 

NRPM 

1 . 0000000 

35 

MAXEP 

155000.0000000 

16 

RPMN 

2600.0000000 

36 

NBLAD 

3.0000000 

17 

DELRPM 

. 0000000 

37 

PDIA 

1 . 9000000 

18 

WEIGHT 

1250.0000000 

38 

WSET 

.0175000 

19 

NXCG 

3.0000000 

39 

TPSET 

- . 0036000 

20 

xcc  X 

25.0988597 

40 

ETAG 

.0000000 

TO  CHANGE  VALUE  TYPE  I , NEW  F 
4  0- 
7  100. 

0  0 

TYPE  -1,0,1  TO  REVIEW, TERMINATE  OR  VIEW  REMAINING  DATA 
-1 


»*“  FTSIM  SETUP  DATA-TO  CHANGE  TYPE  I.F  **** 


-  TERMINATE  CHANGES  WITH  0,0  - 


#(I> 

QUANTITY 

VALUE (F) 

#(I> 

QUANTITY 

VALUE (F) 

1 

TTOT 

. 0000000 

21 

DELXCG 

.1135060 

2 

TSAMP 

. 1000000 

22 

ZCG 

.0591000 

3 

NH 

1.0000000 

23 

PEFF  (1) 

1 . 0000000 

4 

HN 

. 0000000 

24 

PEFF (2) 

1 . 0000000 

5 

DEI, II 

, 0000000 

25 

PEFF  (3) 

1 . 0000000 

6 

NV 

21,0000000 

26 

PEFF  (4) 

1 . 0000000 

7 

VN 

100.0000000 

27 

PEFF  (5) 

1 . 0000000 

8 

DELV 

3 . 0000000 

28 

PEFF (6) 

1 . 0000000 

9 

NPHI 

1 . 0000000 

29 

PEFF (7) 

. 0000000 

10 

PHIN 

. 0000000 

30 

XT 

3 . 0000000 

11 

DELPHI 

31 

ZT 

- . 5500000 

12 

NPL 

. 0000000 

32 

XTH 

-3.6600001 

13 

PLN 

. 0000000 

33 

ZTH 

. 0000000 

14 

DELPL 

. 0000000 

34 

THSET 

. 0000000 

15 

NRPM 

1 . 0000000 

35 

MAXEP 

155000 . 0000000 

16 

RPMN 

2600 . 0000000 

36 

NBLAD 

3 . 0000000 

17 

DELRPM 

. 0000000 

37 

PDIA 

1 . 9000000 

18 

WEIGHT 

1250.0000000 

38 

WSET 

.0175000 

19 

NXCG 

3 . 0000000 

39 

TPSET 

- . 0036000 

20 

XCG  % 

25.0988597 

40 

ETAG 

. 0000000 

TO  CHANCE  VALUE  TYPE  I, NEW  F 

0  0 

TYPE  -1,0, 1  TO  REVIEW, TERMINATE  OR  VIEW  REMAINING  DATA 
1 


“**  FTSIM  SETUP  DATA-TO  CHANGE  TYPE  I,F  **** 


-  TERMINATE  CHANGES  WITH  0,0 


#(I) 

QUANTITY 

VALUE (F) 

#(i> 

QUANTITY 

VALUE (F) 

1 

MASS 

1250.0000000 

28 

A2 

2 . 2000000 

2 

IXX 

1355.0000000 

29 

A3 

. 0000000 

3 

IYY 

2466.0000000 

30 

80 

« 0000000 

4 

IZZ 

3660 . 0000000 

31 

81 

» 0000000 

5 

IXZ 

• 0000000 

32 

82 

- . 4200000 

6 

SW 

15 . 0000000 

33 

83 

- . 1200000 

7 

cw 

1 . 5700001 

34 

CD0T 

•  0000000 

8 

BW 

10  > 0000000 

35 

CDLT 

. 0000000 

9 

ST 

2.7200000 

36 

CMT0 

. 0000000 

10 

CT 

. 8700000 

37 

CMQW 

. 0000000 

11 

CETA 

. 3370000 

38 

CLP 

- . 4050000 

12 

XP 

-1.2660000 

39 

CLXI 

-  .  3600000 

13 

ZP 

. 2600000 

40 

KWING 

. 0000000 

14 

CWP 

2 . 0000000 

41 

KFUSE 

. 0000000 

15 

BFW 

. 8100000 

42 

NPSFI 

. 0100000 

16 

BTAIL 

3.1099999 

43 

NPSFR 

. 0000000 

17 

CL0 

.225610?) 

44 

MPSFI 

. 0000000 

18 

CLAL 

A. 0700002 

45 

MPSFR 

. 0000000 

19 

CD0 

.0380000 

46 

CDB 

.0125000 

20 

CDAL 

.0650000 

47 

GTR 

1 . 0000000 

21 

CM0 

- .0250000 

48 

BL0 

25.0000000 

22 

CMAL 

. 0000000 

49 

WTR 

. 4400000 

23 

EPS0 

. 0000000 

50 

TAPERF 

. 5000000 

24 

EPSAL 

. 4000000 

51 

XWC 

- 1 . 2000000 

25 

QTOQ 

1 . 0000000 

52 

ZWC 

.  2600000 

26 

A0 

*  0000000 

53 

S0S 

1 . 5000000 

27 

A1 

3.1600001 

54 

XQARTC 

- 1 . 2000000 

TO  CHANGE  VALUE  TYPE  I , NEW  F 

0  0 

TYPE  -1.0,1  TO  REVIEW, TERMINATE  OR  VIEW  REMAINING  DATA 

0 


APPENDIX  4.  EXAMPLE  OE  TIME  HISTORIES  ANALYSIS 


SDOFAP. SETUP 


"*»*  ACSL.  COMMAND  FILE  TO  PERFORM  TIME  HISTORY  ANALYSIS.  *»*" 

"***  APPLICATION:  LONGITUDINAL  STABILITY  STUDY  OF  A  LIGHT  ***" 

"***  AIRCRAFT  INFLUENCED  BY  POWER  EFFECTS.  ***" 

S  PRN=9,TCWPRN=72.RRR=21 

"**»  TIME  HISTORY  ANALYSIS  OF  THE  SHORT  PERIOD  RESPONSE  ***" 

"»**  OF  A  LIGHT  AIRCRAFT  TO  AN  ELEVATOR  PULSE  INPUT.  ***" 

S  TSTP=5  .  0,  NSTP=1, CINT=  .05 

S  CALPLT=.F. ,PRNPLT=.F. ,STRPLT=.T. ,CRDSPL=.T. 

OUTPUT  VEK, ALPHAD, Q, ALT, ETAD, GAMMAD, THE TAD , TIME , AN, 'NCIOUT’=10 
PREPAR  TIME, VEK, ALPHAD, Q, ETAD, ALT, GAMMAD, THETAD. AN 

S  DPL=1. 

S  THRMAX=0 . 0 , TSTART=0 . 0, TPULSE=TSTP , TREPET=200 . 

S  ETAMAX= -5.0, ESTART=0 . 01 , EPULSE= , 5 , EREPET=200 . 

S  TST0P=TSTP,DXCCP=7. 22968 
S  DALT0= 10000. 

START 

D  CGPOS, PLS 
S  CMD=DIS 

S  RRR=21 

S  TITLE= '  LIGHT  AIRCRAFT  -  SHORT  PERIOD  MODE. 

RANGE  ETAD 
S  CMD=DIS 

PLOT  ETAD, 'XHI ’=TSTP, 'HI '=HI1, 'L0'=L01, 'XTAG'=' (SEC) ' , ' TAG ' = ' (DEG) ' 

S  TITLE= '  ' 

RANGE  ALPHAD, Q 
S  CMD=DIS 

PLOT  ALPHAD, 'HI '=HI1, ’L0'=L01, ' TAG ' = ' (DEG) '  ... 

,Q. 'HI'=HI2. ' LO ' =L02 , ' TAG ' = ' (RAD/SEC) ' 

RANGE  AN, VEK 
S  CMD=DIS 

PLOT  AN, 'HI ' =HI 1 , 'L0’=L01, ' TAG ' = ' (  G  ) '  ... 

, VEK , 'HI '=HI2, 'L0'=L02, ' TAG ' = ' (KTS) ' 

S  TITLE= '  LIGHT  AIRCRAFT  -  SHORT  PERIOD  MODE. 

RANGE  GAMMAD, THETAD 
S  CMD=DIS 

PLOT  GAMMAD, 'HI '=HI1, 'L0'=L01, ' TAG ' = ' (DEC) '  ... 

, THETAD, 'HI '=HI2. 'L0'=L02, ' TAG ' = ' (DEG) ’ 
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S  CMD=DIS 


TIKE  HISTORY  ANALYSIS  OF  THE  SHORT  PERIOD  RESPONSE  ***' 

OF  A  LIGHT  AIRCRAFT  TO  AN  ELEVATOR  PULSE  INPUT.  ***' 

S  TSTP=5.0,NSTP=1.CINT=.05 

S  CALPLT= . F . . PRNPLT= . F . , STRPLT= . T . . CRDSPL= . T . 

OUTPUT  VEK, ALPHAD, Q, ALT, ETAD, GAMMAD, THETAD, TIME, AN, 'NCIOUT'=10 
PREPAR  TIME, VEK, ALPHAD, Q, ETAD, ALT, GAMMAD, THE TAD. AN 

S  DPL=1. 

S  THRMAX=0 . 0 , TSTART=0 . 0 , TPULSE=TSTP , TREPET=200 . 

S  ETAMAX=- 5 .0, ESTART=0 . 01 . EPULSE= . 5 , EREPET=200 . 

S  TST0P=TSTP,DXCGP=7. 22968 
S  DALT0=10000. 

START 


VEK 
ALT 
THE TAD 

99.9999059 
10000 . 0000 
4.49191191 

ALPHAD  2.36140201 
ETAD  0.80067985 
TIME  0. 

Q  0. 

GAMMAD  2.13050990 
AN  0.99695606 

VEK 
ALT 
THE TAD 

99.7383439 

10004.2341 

11.6790795 

ALPHAD 

ETAD- 

TIME 

8.08575800 

■4.19923639 

0 . 50000000 

Q  0.40193779 
GAMMAD  3.59332151 
AN  1.84753308 

VEK 
ALT 
THE TAD 

98.4150283 

10014.5120 

16.2218415 

ALPHAD 

ETAD 

TIME 

7.82560635 

0.80067985 

1.00000000 

Q-0. 00434133 
GAMMAD  8.39623514 
AN  1.79180178 

VEK 
ALT 
THE  TAD 

97.0084709 

10030.8854 

14 . 9046647 

ALPHAD 

ETAD 

TIME 

4.13979811 

0.80067985 

1 . 50000000 

Q-0. 05634470 
GAMMAD  10.7648666 
AN  1.19691563 

VEK 
ALT 
THE  TAD 

95.7332492 

10049.0153 

13.7025095 

ALPHAD 

ETAD 

TIME 

2.60163220 

0.80067985 

2 • 00000000 

Q-0. 02671067 
GAMMAD  11.1008773 
AN  0.94981021 

VEK 
ALT 
THE  TAD 

94.5503649 

10066.8256 

13.2415065 

ALPHAD 

ETAD 

TIME 

2.43318653 

0.80067985 

2 . 50000000 

Q-0. 00849723 
GAMMAD  10.8083200 
AN  0.90641519 

VEK 
ALT 
THE TAD 

93.4366668 

10083.8828 

13.0790781 

ALPHAD 

ETAD 

TIME 

2.62361519 

0.80067985 

3 . 00000000 

Q-0. 00441509 
GAMMAD  10.4554629 
AN  0.91313888 

VEK 
ALT 
THE  TAD 

92 . 3810440 
10100.2189 
12.9373689 

ALPHAD 

ETAD 

TIME 

2.80481571 

0.80067985 

3 . 5 0000000 

Q-0. 00586753 
GAMMAD  10.1325532 
AN  0.91809176 

VEK 
ALT 
THE TAD 

91.4306631 

10115.8837 

12.7424994 

ALPHAD 

ETAD 

TIME 

2 . 92349396 
0.80067985 

4 . 00000000 

Q-0. 00661579 
GAMMAD  9.81900547 
AN  0.92266517 

VEK 
ALT 
THE TAD 

90.6550760 

10130.9422 

12 . 5656171 

ALPHAD 

ETAD 

TIME 

3.04215094 

0.80067985 

4 . 50000000 

Q-0. 00630872 
GAMMAD  9.52346617 
AN  0.92371201 

VEK 

ALT 

THETAD 

89 . 9308754 
10145  4264 

12 . 3637042 

ALPHAD 

ETAD 

TIME 

3 . 13963864 
0.80067985 

5 . 00000000 

Q-0 .00790029 
CAMMAD  9.22406554 
AN  0.92242613 

PLS  1.00000000 


D  CGPOS.PLS 

CGPOS  32.3285397 
S  CMD=DIS 
S  CMD=10 

S  RRR=21 

S  TITLE='  LIGHT  AIRCRAFT  -  SHORT  PERIOD  MODE.  ' 

RANGE  ETAD 

ETAD-4. 19923639  0.80067985 

S  CMD=DIS 
S  HI 1=2 . , L01=-6. 

S  CMD=10 

PLOT  ETAD, 'XHI’=TSTP, 'HI'=HI1, 'LO'=L01, ' XTAG ' = ' (SEC) ' TAG ' = ' (DEG) ' 

S  TITLE=*  ' 

RANGE  ALPHAD.Q 

ALPHAD  2.36140201  9.38391014 

Q-0. 05862956  0.40193779 

S  CMD=DIS 

S  HI 1=14. ,L01=-2. ,HI2=.6,L02=- . 2 
S  CMD=10 

PLOT  ALPHAD,  'HI’=HI1,  'L0'=L01,  'TAG'='  (DEG)  *  ... 

,Q, 'HI ' =HI 2 , ' L0'=L02, ' TAG ' = 1 (RAD/SEC) ' 

RANGE  AN, VEK 

AN  0.90631505  2.08298743 

VEK  89.8613617  100.000728 

S  CMD=DIS 

S  HI 1=2 . 5 , L01= . 5, HI 2=105 . , L02=85 . 

S  CMD=10 

PLOT  AN, ' HI ' =HI 1 , ' LO ' =L01 , ' TAG ' = ' (  G  ) ’  ... 

.VEK, 1 HI ' =HI 2 , 1 L0'=L02, ' TAG 1 = ' (KTS) ' 

S  TITLE= '  LIGHT  AIRCRAFT  -  SHORT  PERIOD  MODE.  ' 

RANGE  GAMMAD , THE TAD 

GAMMAD  2.09020455  11.1164848 

THETAD  4.49191191  16.2218415 

S  CMD=DIS 

S  HI 1=12 . , L01=2 . , HI 2=20. , LO2=0 . 

S  CMD=10 

PLOT  GAMMAD. 'HI'=HI1. 'L0'=L01, ' TAG ' = ' (DEG) '  ... 

.THETAD. 'HI’=HI 2. 'L0'=L02, ' TAG ' = ' (DEG) ' 

S  CMD=DIS 
STOP 


* 


■* 


*.*-•.  .  ■ 


SDOFAP . SETUP 


*****  ACSL.  COMMAND  FILE  TO  PERFORM  TIME  HISTORY  ANALYSIS.  ***" 

"***  APPLICATION:  LONGITUDINAL  STABILITY  STUDY  OF  A  LIGHT  ***** 

"***  AIRCRAFT  INFLUENCED  BY  POWER  EFFECTS.  ***** 

S  PRN=9,TCWPRN=72,RRR=21 

"**•  TIME  HISTORY  ANALYSIS  OF  THE  PHUGOID  RESPONSE  OF  ***" 

*****  A  LIGHT  AIRCRAFT  TO  AN  ELEVATOR  PULSE  INPUT.  *»*" 

S  TSTP=120.p.NSTP=l,CINT=.05 
S  CALPLT= . T . , PRNPLT= . F . , STRPLT= . F . , GRDCPL= . T . 

OUTPUT  VEK, ALPHAD, Q, ALT, ETAD, GAMMAD, THE TAD, TIME , AN , *NCIOUT'=240 
PREPAR  TIME, VEK, ALPHAD, Q, ETAD, ALT, GAMMAD, THETAD, AN 

S  DPL=0. 

S  THRMAX=p .0, TSTART=0 .0, TPULSE=TSTP , TREPET=200 . 

S  ETAMAX=1  .0,  ESTART=0.0.l ,  EPULSE=5  .  ,  EREPET=200 . 

S  TSTOP=TSTP , DXCGP=7 . 22968 
S  DALT0=10000. 

START 

D  CGPOS , PLS 
S  CMD=DIS 

S  RRR=21 

S  TITLE= '  LIGHT  AIRCRAFT  -  PHUGOID  MODE.  ’ 

RANGE  ETAD 
S  CMD=DIS 

PLOT  ETAD. ’ XHI ' =TSTP , 'HI*=HI1, 'LO'=LOI, 'XTAG'=' (SEC) ’ TAG ’ = ' (DEG) ' 

S  TITLE='  ' 

RANGE  ALPHAD. Q. AN 
S  CMD=DIS 

PLOT  ALPHAD. *HI'=HU,  *LO'=LOI. 'TAG'=* (DEG) *,Q,  *HI’=HI2. 'L0'=L02  . . . 

, 'TAG’=* (RAD/SEC) ’, AN, 'HI '=HI3, 'L0'=L03, *TAG'=' (  G  )' 

RANGE  VEK. GAMMAD, THETAD 
S  CMD=DIS 

PLOT  VEK, 'HI '=HI1, 'L0'=L01, *TAC'=' (KTS) ’ .GAMMAD, 'HI *=HI2, 'L0’=L02  . 
, ' TAG ' = ' (DEG) '.THETAD, *HI'=HI3, ’L0'=L03, *TAG'=*  (DEG) ’ 


TIME  HISTORY  ANALYSIS  OF  THE  PHUGOID  RESPONSE  OF  ***' 

' *  *  *  A  LIGHT  AIRCRAFT  TO  AN  ELEVATOR  PULSE  INPUT.  ***' 

S  TSTP=120.0,NSTP=1,CINT=.05 
S  CALPLT= . T . . PRNPLT= . F . , STRPLT= . F . , GRDCPL= . T . 

OUTPUT  VEK, ALPHAD, Q, ALT, ETAD, GAMMAD, THETAD, TIME , AN. ’NCIOUT'=240 
PREPAR  TIME, VEK, ALPHAD, Q, ETAD, ALT, GAMMAD, THETAD, AN 

S  DPL=0. 

S  THRMAX=0 .  0 , TSTART=0 . 0 , TPULSE=TSTP , TREPET=200 . 

S  ETAMAX=1.0,ESTART=0.01,EPULSE=5.  ,EREPET=200. 

S  TSTOP=TSTP, DXCGP=7 . 22968 
S  DAL  T0= 10000. 

START 


VEK  99.9999059 

ALT  1 0000.0000 
THETAD-3. 43137409 

ALPHAD  2.69495477 
ETAD  0.49637286 
TIME  0. 

Q  0. 

GAMMAD-6. 12632886 
AN  0.99929857 

VEK  121.937459 

ALT  9399.40124 
THETAD-3. 45999241 

ALPHAD 

ETAD 

TIME 

1.93136242 

0.49637286 

12 .0000000 

Q  0.04324758 
GAMMAD- 5. 39135482 
AN  1.31508105 

VEK  87.5520883 

ALT  9521.73771 
THETAD  6.23735405 

ALPHAD  3.43412635 
ETAD  0.49637286 
TIME  24.0000000 

Q-0. 02549850 
GAMMAD  2.80322770 
AN  0.85162012 

VEK  99.1196767 

ALT  9218.32978 
THETAD- 12. 6037528 

ALPHAD  2.65628993 
ETAD  0.49637286 
TIME  36.0000000 

Q-0. 00102716 
GAMMAD- 15 . 2600428 
AN  0.97578389 

VEK  109.290973 

ALT  8795.32924 
THETAD  0.02009332 

ALPHAD  2.31813959 
ETAD  0.49637286 
TIME  48.0000000 

Q  0.01964675 
GAMMAD- 2. 29804627 
AN  1.12601581 

VEK  91.1228162 

ALT  8763.29354 
THETAD- 1.04949672 

ALPHAD 

ETAD 

TIME 

3.19802476 

0.49637286 

60 . 0000000 

Q-0. 01958144 
GAMMAD-4. 24752148 
AN  0.89270209 

VEK  103.209102 
ALT  8414.01146 
THETAD-7. 7325547 2 

ALPHAD 

ETAD 

TIME 

2 .51945226 
0.49637286 
72.0000000 

Q  0.00702534 
GAMMAD- 10 . 2520070 
AN  1.03622142 

VEK  102.264847 
ALT  8150.97446 
THETAD-0. 12452573 

ALPHAD 

ETAD 

TIME 

2.59925787 

0.49637286 

84 . 0000000 

Q  0.00533936 
GAMMAD-2 . 72378360 
AN  1.03013203 

VEK  95.5421614 
ALT  8000.27666 
THETAD- 3. 90601896 

ALPHAD 

ETAD 

TIME 

2 . 92362864 
0.49637286 

96 . 0000000 

Q-0. 00980669 
GAMMAD-6. 82964760 
AN  0.94358945 

VEK  103.120145 

ALT  7681.26595 
THETAD-4 .74415664 

ALPHAD 

ETAD 

TIME 

2 . 54128177 

0 . 49637286 

108 . 000000 

Q  0.00674860 
GAMMAD -7 . 28543841 
AN  1.03801802 

VEK  99.3821384 
ALT  7470.98509 
THETAD- 1 .43838404 

ALPHAD 

ETAD 

TIME 

2 . 73222594 

0 . 49637286 

1 20 . 000000 

Q-0. 00101423 
GAMMAD-4 . 17060998 
AN  0.99260128 

D  CCPOS.PLS 

CGPOS  32.3285397  PLS  0. 

S  CMD=DIS 
S  CMD=10 

S  RRR=21 

S  TITLE='  LIGHT  AIRCRAFT  -  PHUGOID  MODE.  ' 

RANGE  ETAD 

ETAD  0.49637286  1.49635611 

S  CMD=DIS 
S  HI 1=2 . , LO1=0 . 

S  CMD=10 

PLOT  ETAD. 'XHI'=TSTP, 'HI'=HI1, 'LO'=L01, 'XTAG'  =  ' (SEC)  *  TAG ' =  (DEG) 

S  TITLE= ' 

RANGE  ALPHAD, Q, AN 

ALPHAD-0. 23557132  3.77131611 

Q-0. 08034203  0.05005342 

AN  0.61702932  1.32090147 

S  CMD=DIS 

S  HI 1=4 . , L01=-l . , HI2= . 1, L02=- . l.HI3=l . 4, L03= . 4 
S  CMD=10 

PLOT  ALPHAD.  ’HI'=HI1.  'L0'=L01,  ' TAG '  = '  (DEG)  \Q,  ,HI'=HI2.  'L0'=L02  .  .  . 

, ' TAG '  = '  (RAD/SEC)  ' . AN . ' HI  1 =HI 3 , '  LO ' =L03 ,  ' TAG  *  =  * (  G  )' 

RANGE  VEK . GAMMAD , THE TAD 

VEK  03.2781743  122.456445 

CAMMAD- 20. 0404889  6.21916033 

THETAD-19. 6957197  8.95900022 

S  CMD=DIS 

S  HI1=130. ,LO1=80. .HI2=15. ,L02=-35. , HI 3=25.  ,L03=*25. 

S  CMD=10 

PLOT  VEK. 'HI ' =HI1. 'L0’=L01, ' TAG ' = ' (KTS) '.GAMMAD, 'HI’=HI2, 'L0’=L02  . . . 
,  ’ TAG 1  =  ' (DEG)  '.THETAD,  'HI'=HI3,  ’L0'=L03,  ' TAG '  = ' (DEG)  ' 

S  CMD=DI S 
STOP 


TIME  .jfl2  (SEC) 


% 


SDOFAP .  SETUP 


"***  ACSL.  COMMAND  FILE  TO  PERFORM  EIGEN  ANALYSIS.  ***" 
"***  APPLICATION:  LONGITUDINAL  STABILITY  STUDY  OF  ***" 
"***  A  LIGHT  AIRCRAFT  INFLUENCED  BY  POWER  EFFECTS. 


S  PRN=9,TCWPRN=72 

"***  EIGEN  ANALYSIS  OF  LONGITUDINAL  MOTION  OF  A  ***" 
"***  LIGHT  AIRCRAFT  INFLUENCED  BY  POWER  EFFECTS.  *»*" 

S  NSTP=1,CINT=.05 
S  DXCCP=7. 22968, TSTOP=p. 

S  DALT0=lf)000 . 

PREPAR  TIME , VEK , ALPHAD, Q, ETAD , ALT , GAMMAD , THE TAD, AN 


”***  EIGEN  ANALYSIS  ***" 

S  DPL=p. 

START 

D  CGPOS.PLS 

ANALYZ  ' FREEZE ' =X, Y. Z. P , R , BETAR , TAU1 ,  TAU3 
ANALYZ  'EIGVEC'=.T. , ’EIGEN' 


"***  EIGEN  ANALYSIS  «**" 

S  BEGIN= . F . 

S  DPL=I . 

START 

D  CGPOS , PLS 


ANALYZ  'EIGVEC'=.T.  ’EIGEN 


SDOFAP.L 


'***  EIGEN  ANALYSIS  OF  LONGITUDINAL  MOTION  OF  A  ***  ' 
’***  LIGHT  AIRCRAFT  INFLUENCED  BY  POWER  EEFECTS .  ***' 

S  NSTP=1,CINT=.05 
S  DXCGP=7. 22968,  TSTOP=0 . 

S  DALT0=10000. 

PREPAR  TIME , VEK , ALPHAD, Q, ETAD, ALT, GAMMAD , THE TAD, AN 


'***  EIGEN  ANALYSIS  ***' 

S  DPL=0. 

START 

D  CGPOS , PLS 

CGPOS  32.3285397  PLS  0. 

ANALYZ  *  FREEZE ' =X , Y , Z , P . R . BETAR , TAU1 , TAU3 
ANALYZ  'EIGVEC'=.T. , 'EIGEN' 


COMPLEX  EIGEN  VALUES  IN  ASCENDING  ORDER 

REAL  IMAGINARY  FREQUENCY 

1  -3.2785E-16 

2  -0.01886002  18552241  0.186479 

4  -1.84179949  ♦/- 1 . 71229370  2.514791 


DAMPING 

0.101138 

0.732387 


COMPLEX  EIGEN  VECTORS 
1 

1  -0.9819847  0. 

2  0.0290791  0. 

3  1 . 461E-16  0. 

4  0.1867097  0. 

5  2 . 800E-07  0. 


2 

2 . 849E-04- 3 . 806E-05 
0.0095129-0.0012707 
1.127E-04  0.0035793 
0.0142453  0.9998451 
7. 715E-05-0. 0013708 


3 

2 . 849E-04  3.806E-05 
0.0095129  0.0012707 
1 . 127E-04-0. 0035793 
0.0142453-0.9998451 
7 . 715E-05  0.0013708 


4 

0.0025064  0.0031429 
0.0836763  0.1049257 
0.0511200-0.6733628 
0.4396764  0.4369320 
0.3589739  0.1133794 


5 

0.0025064-0.0031429 
0.0836763-0.1049257 
0.0511200  0.6733628 
0.4396764-0.4369320 
0.3589739-0.1133794 


'***  EIGEN  ANALYSIS  ***' 


S  BEGIN= . F . 
S  DPL=1. 


START 

D  CGPOS. PLS 

CGPOS  32.3285397  PLS  1 .00000000 


ANALYZ  'EICVEC'=.T.  'EIGEN' 


COMPLEX  EIGEN  VALUES  IN  ASCENDING  ORDER 

REAL  IMAGINARY  FREQUENCY 

1  -6  9337E-17 

2  -0.01700625  V-0. 15272301  0.153667 

4  -1.99230691  +/-1 .74483906  2.648349 


DAMPING 

0.110670 

0.752283 


COMPLEX  EIGEN  VECTORS 
1 

1  -0.9992318  0. 

2  -0.0391893  0. 

3  -1.728E-16  0. 

4  -4 . 805E-09  0. 

5  4.605E-12  0. 


2 

-1 . 808E-04  2 . 549E-04 
0 . 0046099-0 . 0064994 
0.0018298  0.0016304 
0.6663900  0.7455544 
-0.0012826-0.0013281 


3 

-1 .808E-04-2 . 549E-04 
0.0046099  0.0064994 
0.0018298-0.0016304 
0 . 6663900-0 . 7455544 
-0.0012826  0.0013281 


4 

3. 778E-04-0. 0052915 
-0.0096333  0.1349206 
0 . 5096059-0 . 5043772 
-0.0298649  0.5624622 
0.1840649  0.3412410 


5 

3 . 778E-04  0.0052915 
-0.0096333-0.1349206 
0.5096059  0.5043772 
-0.0298649-0.5624622 
0.1840649-0.3412410 


STOP 


2  WORDS  TABLE  SPACE  USED 


APPENDIX  6.  EXAMPLE  OF  JACOBIAN  ANALYSIS 


SDOFAP.SE TUP 


"**»  ACSL.  COMMAND  FILE  TO  PERFORM  JACOBIAN  ANALYSIS. 
"***  APPLICATION:  LONGITUDINAL  STABILITY  STUDY  OF  A 
"***  LIGHT  AIRCRAFT  INFLUENCED  BY  POWER  EFFECTS. 


S  PRN=9 , TCWPRN=7  2 

"***  JACOBIAN  ANALYSIS  TO  DEDUCE  NON-DIMENSIONAL  ***" 
"***  AERODYNAMIC  DERIVATIVES  OF  A  LIGHT  AIRCRAFT 
"***  INFLUENCED  BY  POWER  EFFECTS.  ***” 

S  NSTP=1,CINT=.05 
S  DXCGP=7 .22968, TSTOP=0 . 

S  DAL  10=1000(3. 

PREPAR  TIME , VEK, ALPHAD, Q, ETAD, ALT, GAMMAD , THE TAD, AN 


"***  JACOBIAN  ANALYSIS  ***" 

S  DPL=1. 

START 

D  CGPOS.PLS 

ANALYZ  ' FREEZE ' =X, Y. Z , P . R , BETAR , TAU1 , TAU3 
ANALYZ  'JACOB' 

S  PRN=DIS 

ANALYZ  'EIGEN'  $"***  THIS  COMMAND  IS  USED  ONLY  TO  * 
$"***  INITIATE  THE  JACOBIAN  ANALYSIS  * 


S  PRN=9 

"»*«  JACOBIAN  ANALYSIS  ***" 

S  BEGIN=.F. 

S  DPL=1. 

START 

D  CGPOS.PLS 
ANALYZ  'JACOB' 


***" 
***» 
*  *  »  « 


*« 
*  » 


S  PRN=DIS 
ANALYZ  'EIGEN' 


SDOFAP.L 


'***  JACOBIAN  ANALYSIS  TO  DEDUCE  NON-DIMENSIONAL  ***' 
'***  AERODYNAMIC  DERIVATIVES  OF  A  LIGHT  AIRCRAFT  ***' 
'***  INFLUENCED  BY  POWER  EFFECTS. 

S  NSTP=1,CINT=.05 
S  DXCCP=7.22968,TSTOP=0. 

S  DALT0=10000. 

PREPAR  TIME , VEK , ALPHAD , Q , ETAD , ALT , GAMMAD , THE TAD , AN 


’***  JACOBIAN  ANALYSIS  ***' 


S  DPL=1 . 


START 

D  CGPOS.PLS 

CGPOS  32.3285397 


PLS  1 .00000000 


ANALYZ  ' FREEZE ' =X. Y. Z , P , R, BETAR , TAU1 , TAU3 
ANALYZ  'JACOB' 


ROW  VECTOR  NAMES 

TAU0  1 

VT  4 

COLUMN  VECTOR  NAMES 
TAU0DT  1 

DVT  4 


TAU2 

ALPHAR 


TAU2DT 

DALPHR 


MATRIX  ELEMENTS  -  ROWS  ACROSS.  COLUMNS  DOWN 


1  2  3  4  5 

1  0.  0.  -0.0195946  0.  0. 

2  0.  0.  0.4996159  0.  0. 

3  -3 . 991E-04  0.0101762  -2.5035555  -9.349E-05  -3.3652607 

4  0.7681246  -19.585314  -0.0174290  -0.0431410  4.6697616 

5  4.730E-04  -0.0120811  0.9798915  -0.0051029  -1.4719298 


S  PRN=DIS 


***  DIMENSIONAL  JACOBIAN. 


*\  • 

Jii . 

A', '•'a' 
A'aV'A' 


TO 


1  v  V  V  v 
• 

S  V  «.\s' 

a’aIn  -s  . 


■V-.-VV 

A  A  t*  .  j*  j 


****''-4mm 

--a' A".'-'. 

'-A.-v.; 

*  '  .  ■  ,  *  A* 

> _ i 

A>y-> 

V  .' 


>  a  a 

V  \  . 
--  .*  .* 

•  .  •  A  A  A 
-  .  -'A  A  Vi 
.'A'A'.'.'A 

•a'.Va'j 


•••r.-v-y 


-.0431410  4.6697616 

-.0051029  -1.4719298 


.  0000000 


-.0174290  -19.5853136 


.  9798915 


-.0000935  -3.3652607  -2.5035555 


.4996159 


.0120811 
.0101762 
.  00000^30 


NOTE:  ASSUMED  RELATIONSHIPS  IN  USE  FOR  CALCULATION 
OF  ANGULAR  RATE  DERIVATIVES. 


*  A  ..  A  A 


“N'aVV 

•>A>. 

AaVV 

vA>YV 

■>'a;.v 

kJU  •»*.’« 


"."a::;-: 


M  Q 


***  NON-DIMENSIONAL  JACOBIAN.  *** 


-  .0005653 
- . 0040057 
- . 0^00010 
. 0000000 


.0010215 
- .0192886 
- .0005779 
. 0000000 


- .0174290 
. 9798915 
-  .0328073 
.4996159 


-.0021422 
- .0000792 
. PPPPP09 
*  PPPP^PP 


NON-DIMENSIONAL  AERODYNAMIC  DERIVATIVES.  *** 


CTV 

CLV 

CDV  = 

CMV  = 

CLALPHA  = 
CDALPHA  = 
CMALPHA  = 
CLDALPHA  = 
CMDALPHA  = 
CLQ  = 

CMQ 


- .2208641 
- .0606333 
- .0509233 
- .0150985 
4.4965035 
.2632860 
- .2876742 
1.3753121 
-3.6791788 
3.2512172 
-8.6975235 


JACOBIAN  ANALYSIS  ***' 

S  BEGIN= . F . 

S  DPL=1. 

TART 

CGPOS  PLS 

CGPOS  32.3285397  PLS  1 . 00000000 


ANALYZ  ' JACOB 


ROW  VECTOR  NAMES 

TAU0  1  TAU2  2  Q  3 

VT  4  ALPHAR  5 

COLUMN  VECTOR  NAMES 

TAU0DT  1  TAU2DT  2  QDOT  3 

DVT  4  DALPHR  5 

MATRIX  ELEMENTS  -  ROWS  ACROSS,  COLUMNS  DOWN 

12  3  4 

1  {3.  0.  -(3.0195946  0. 

2  0.  0.  0.4996159  0. 

3  -3 . 991E-04  0.0101762  -2.5035555  -9.349E-05 

4  0.7681246  -19.585314  -0.0174290  -0.0431410 

5  4 . 738E-04  -0.0120811  0.9798915  -0.0051029 


S  PRN=DIS 


**«  DIMENSIONAL  JACOBIAN.  *** 


-.0431410  4.6697616  -.0174290  -19.5853136 

-.0051029  -1.4719298  .9798915  -.0120811 

-.0000935  -3.3652607  -2.5035555  .0101762 

.0000000  .0000000  .4996159  .0000000 


NOTE:  ANGULAR  RATE  DERIVATIVES  CALCULATED  FROM 
ELEMENTS  OF  THE  JACOBIAN  MATRIX. 


BEWARE  OF  POSSIBLE  INACCURACY  ASSOCIATED  WITH  SMALL  CLIMB  ANGLES. 


*•*  NON-DIMENSIONAL  JACOBIAN.  **» 


-.0005653  .0010215  -.0174290  -.0021422 

-.0040057  -.0192886  .9798915  -.0000792 

-.0000010  -.0005779  -.0328073  .0000009 

.0000000  .0000000  .4996159  .0000000 


***  NON-DIMENSIONAL  AERODYNAMIC  DERIVATIVES 


CTV  = 

CLV 

CDV 

CMV 

CLALPHA  = 
CDALPHA  = 
CMALFHA  = 
CLDALPHA  = 
CMD ALPHA  = 
CLQ 

CMO  = 


-.22(58641 
- .0616025 
-.0509233 
- .0169414 
4.4918363 
.2632860 
- .2965481 
1.1333446 
-4.1392378 
3.6067239 
-8.2467156 


SUBROUTINE  INTERM (A) 


»****• 
****** 
*****  *  * 


PROVIDES  JACOBIAN  REDUCTION  ROUTINE  * 
FLIGHT  CONDITION  DATA. 
THIS  SUBROUTINE  MUST  BE  APPENDED  TO 


ZZREDC '  WITH 
SDOFAP . ACSL 


******* 

******* 

******* 


INTEGER  ROCOL (5) 

DIMENSION  A (5, 5) ,AA(5, 5) .AAA (5, 5) 
DIMENSION  B (15) 

LOGICAL  LINEAR 
CHARACTER* 1  ANS 

DATA  ROCOL/4,5,3,2,1/ 

LINEAR  =  .FALSE. 


"***  REARRANGE  ROWS  OF  JACOBIAN.  ***" 

DO  1(3  1=1.5 
DO  2(3  J=1 , 5 
AAA (I , J) =A (ROCOL (I) ,J) 

CONTINUE 

CONTINUE 


"**»  REARRANGE  COLS  OF  JACOBIAN.  ***" 

DO  3(3  1=1,5 
DO  4(3  J=1 , 5 

AA  (I ,  J)  =AAA  (I ,  ROCOL  (J) ) 

CONTINUE 

CONTINUE 


*****  WRITE  DIMENSIONAL  JACOBIAN  TO  FILE.  ***" 
WRITE  (9,8(3) 

DO  70  1=1,4 

WRITE  (9,85)  (AA(I.J) ,J=1,4) 

CONTINUE 


ASCERTAIN  WHETHER  LINEAR  ANALYSIS  IS  TO  BE  USED  ***" 

WRITE  (6.130) 

READ  (5,100)  ANS 

IF  (  ANS  . EQ.  1HY  .OR.  ANS  . EQ .  lHy  ,  LINEAR  =  .TRUE. 

IF  (  LINEAR  )  THEN 
WRITE  (9. 123) 

ELSE 

WRITE (9, 124) 

WRITE (9. 125) 

WRITE (9, 126) 

ENDIF 


80  FORMAT (////, 31H  ***  DIMENSIONAL  JACOBIAN.  ***,//) 

85  FORMAT (/,4X,F11.7,2X,F11.7,2X,F11.7,2X,F11.7) 

100  FORMAT  (Al) 

123  FORMAT (//, IX, 51HNOTE :  ASSUMED  RELATIONSHIPS  IN  USE  FOR  CALCULATION 

, /, IX, 35H  OF  ANGULAR  RATE  DERIVATIVES.  .//) 

124  FORMAT (//.IX, 47HN0TE :  ANCULAR  RATE  DERIVATIVES  CALCULATED  FROM 

,/.lX,38H  ELEMENTS  OF  THE  JACOBIAN  MATRIX.) 

125  FORMAT (//, IX, 4SHBEWARE  OF  POSSIBLE  INACCURACY  ASSOCIATED  WITH 

, 20H  SMALL  CLIMB  ANGLES.  ) 

126  FORMAT  (  1X.45H . - . . . 

/  20H . -  ,//) 

130  FORMAT (//, IX, 46H  SHOULD  THE  ANGULAR  RATE  DERIVATIVES  BE  FOUND  ,/, 
1X.42H  BY  USINC  ASSUMED  RELATIONSHIPS.  (Y/N)  :  ?,4X$) 


200  RETURN 


REDUCE.  F 


SUBROUTINE  ZZREDC (AA, B, LINEAR) 

C  ***  REDUCES  DIMENSIONAL  JACOBIAN  TO  NON  DIMENSIONAL  DERIVATIVES 

PARAMETER  (ISIZE=5, ISIZEM1=4)  !***  SIZE  OF  THE  PROBLEM 

LOGICAL  FLAG,  LINEAR 
CHARACTER *10  C(ll) 

DIMENSION  AA(ISIZE,  ISIZE)  ,B{24)  ,D(H) 

REAL  IYYND.  LT 

DATA  C  /'CTV  = ' , 'CLV  =','CDV  =’.'CMV  =', 

1  'CLALPHA  =  ' , ' CDALPHA  =  ' , ' CMALPHA  = • . • CLDALPHA  = ' , 

1  ' CMDALPHA  = ' , ' CLQ  =  \'CMQ  ='/ 


C  »**  WRITE  NON-DIMENSIONAL  JACOBIAN  TO  FILE  *** 

WRITE  (9, 50) 

10  DO  20  1=1 , ISIZEM1 

WRITE  (9,70)  (AA(I.J) . J=1 , ISIZEM1) 

20  CONTINUE 


C  ***  RECONSTITUTE  FLIGHT  DATA  *** 


SGA 

=B(1) 

CGA 

=B(2) 

SAL 

=B(3) 

CAL 

=B(4) 

CWE 

=B(5) 

CLE 

=B(6) 

CTE 

=B(7) 

CDE 

=B(8) 

AMU 

=B(9) 

IYYND 

=B(10) 

XLT 

=B(H) 

XEPSAL=B (12) 

CW 

=B(13) 

IF  (  LINEAR  )  THEN 

C  ***  APPROXIMATE  SURPLUS  DERIVATIVE  CTV  •** 


□  (1)  =  -  3 . 0*CTE  ! CTV 

C  ***  CALCULATE  NON-DIMENSIONAL  DERIVATIVES  *** 

D (6)  =  CLE-2 .0*AMU‘AA(1, 2)  ! CDALPHA 

D (3)  =  D{1) *CAL*2.0*CWE‘SGA-2  0*AMU*AA{1, 1)  !CDV 

D (11)  =  IYYND* AA ( 3 , 3) / ( 1 . 0*XEPSAL  * AA (2,3))  !CMQ 

D  (9)  =  D(ll) *XEPSA*  ! CMDALPHR 

D (10)  =  -D(ll) *CW/XLT  \ CLQ 

D (8)  =  XEPSAL*D(10)  'CLDALPHA 


*  *  * 

*  *  * 


DENOM  =  2 . 0*AMU*D (8) 


!  DENOMINATOR 


-AA(2, 2) ‘DENOM-CDE  !CLALPHA 

-D  (1)  *SAL-AA(2, 1)  ‘DENOM- 2  .0*CWE*CGA  !CLV 
IYYND*AA(3. 1) *D (9) *AA (2 , 1)  !CMV 

IYYND»AA(3.2) -D (9) *AA(2.2)  ICMALPHA 


“*  APPROXIMATE  SURPLUS  DERIVATIVE  CTV  “* 


D  (1)  =-3.‘CTE 


t  CTV 


CALCULATE  NON-DIMENSIONAL  DERIVATIVES  *“ 


D (8)  =  -CWE*SGA/AA(2,4) -2. ‘AMU  !  CLDA 

DENOM=  2.0‘AMU*D(8)  !  DENC 

D (9)  =  -AA(3,4) *IYYND*DENOM/CWE/SGA  !  CMDA 

D (6)  =  CLE- 2 . ‘AMU*AA (1 , 2)  !  CDAL 

D (5)  =  -AA(2, 2) ‘DENOM-CDE  !  CLAL 

D (10) =  2 . *AMU-AA(2 , 3) ‘DENOM  !  CLQ 

D (3)  =  D(l) *CAL*2.*CWE*SGA-2.*AMU‘AA(1,1)  !  CDV 

D (2)  =  -D (1) ‘SAL -2 . ‘CWE*CGA-DENOM*AA(2 , 1)  !  CLV 

FACT1=  D (9) * (D (1) ‘SAL+D (2) *2 . *CWE‘CGA) /DENOM 

D (4)  =  IYYND*AA(3, 1) +FACT1  !  CMV 

D (7)  =  IYYND*AA(3, 2) +D  (9) * (D (5) +CDE) /DENOM  !  CMALPHA 

D (11) =  IYYND*AA(3, 3)  -D  (9)  *  (2  .  *AMU-D  (10)  )  /DENOM!  CMQ 


!  CLDALPHA 

!  DENOMINATOR 

!  CMD ALPHA 
!  CDALPHA 
!  CLALPHA 
!  CLQ 
!  CDV 
!  CLV 


ENDIF 


“*  OUTPUT  COEFFICIENTS  TO  FILE  SDOFAP.L  “* 

WRITE (9,60) 

DO  30  11=1,11 

WRITE (9. 40)  C(II)  , D  (II) 

CONTINUE 
WRITE (9,80) 

FORMAT (/,4X,A10, 2X.F14.7) 

FORMAT (////, 35H  “*  NON-DIMENSIONAL  JACOBIAN.  “*,//) 

FORMAT (//// 

.  ,50H  “*  NON-DIMENSIONAL  AERODYNAMIC  DERIVATIVES.  »“,//) 

FORMAT (/.4X.F11 .7, 2X , Fll . 7 , 2X , Fll . 7 , 2X , Fll . 7) 

FORMAT (///) 


RETURN 

END 


SUBROUTINE  ACINFO (CTP . SWING, CWINC, ALPT, XLT, XEPSAL) 


INCLUDE  ' FTPAR . F ' 

CTP=CTPROP 
SWI NG=SW 
CWI NG=CW 
ALPT=THSET 


***  CALCULATE  THE  TAILPLANE'S  MOMENT  ARM  (  XLT  ) 
XLT  =  XT  -  XQARTC 


in 


***  THE  RATE  OF  CHANGE  OF  DOWNWASH  WITH  ALPHA  IS  PREFIXED  WITH  *** 
***  AN  'X'  TO  OVERCOME  PROBLEMS  WITH  THE  COMMON  STATEMENT.  '** 

XEPSAL  =  EPSAL 

RETURN 
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